diff --git a/MapControl/Shared/AzimuthalEquidistantProjection.cs b/MapControl/Shared/AzimuthalEquidistantProjection.cs index 1649b54b..1c98c1cc 100644 --- a/MapControl/Shared/AzimuthalEquidistantProjection.cs +++ b/MapControl/Shared/AzimuthalEquidistantProjection.cs @@ -63,13 +63,13 @@ namespace MapControl var c = Math.Acos(p.CosC); var k = c / Math.Sin(c); // p.195 (25-2) - return new Point(EarthRadius * k * p.X, EarthRadius * k * p.Y); // p.195 (22-4/5) + return new Point(EquatorialRadius * k * p.X, EquatorialRadius * k * p.Y); // p.195 (22-4/5) } public override Location MapToLocation(double x, double y) { var rho = Math.Sqrt(x * x + y * y); - var c = rho / EarthRadius; // p.196 (25-15) + var c = rho / EquatorialRadius; // p.196 (25-15) return GetLocation(x, y, rho, Math.Sin(c)); } diff --git a/MapControl/Shared/AzimuthalProjection.cs b/MapControl/Shared/AzimuthalProjection.cs index dc9a3bf1..0979fdc6 100644 --- a/MapControl/Shared/AzimuthalProjection.cs +++ b/MapControl/Shared/AzimuthalProjection.cs @@ -5,9 +5,6 @@ using System.Windows.Media; namespace MapControl { - /// - /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395), p.141. - /// public abstract class AzimuthalProjection : MapProjection { protected AzimuthalProjection() @@ -16,8 +13,6 @@ namespace MapControl Type = MapProjectionType.Azimuthal; } - public double EarthRadius { get; set; } = Wgs84MeanRadius; - public readonly struct ProjectedPoint { public double X { get; } diff --git a/MapControl/Shared/EquirectangularProjection.cs b/MapControl/Shared/EquirectangularProjection.cs index b2832c3b..9644112b 100644 --- a/MapControl/Shared/EquirectangularProjection.cs +++ b/MapControl/Shared/EquirectangularProjection.cs @@ -35,12 +35,12 @@ namespace MapControl public override Point? LocationToMap(double latitude, double longitude) { - return new Point(Wgs84MeterPerDegree * longitude, Wgs84MeterPerDegree * latitude); + return new Point(MeterPerDegree * longitude, MeterPerDegree * latitude); } public override Location MapToLocation(double x, double y) { - return new Location(y / Wgs84MeterPerDegree, x / Wgs84MeterPerDegree); + return new Location(y / MeterPerDegree, x / MeterPerDegree); } } } diff --git a/MapControl/Shared/GnomonicProjection.cs b/MapControl/Shared/GnomonicProjection.cs index f20d3e5f..f03c15c0 100644 --- a/MapControl/Shared/GnomonicProjection.cs +++ b/MapControl/Shared/GnomonicProjection.cs @@ -46,13 +46,13 @@ namespace MapControl var k = 1d / p.CosC; // p.165 (22-3) - return new Point(EarthRadius * k * p.X, EarthRadius * k * p.Y); // p.165 (22-4/5) + return new Point(EquatorialRadius * k * p.X, EquatorialRadius * k * p.Y); // p.165 (22-4/5) } public override Location MapToLocation(double x, double y) { var rho = Math.Sqrt(x * x + y * y); - var c = Math.Atan(rho / EarthRadius); // p.167 (22-16) + var c = Math.Atan(rho / EquatorialRadius); // p.167 (22-16) return GetLocation(x, y, rho, Math.Sin(c)); } diff --git a/MapControl/Shared/Location.cs b/MapControl/Shared/Location.cs index b09c8327..e3ecb9ae 100644 --- a/MapControl/Shared/Location.cs +++ b/MapControl/Shared/Location.cs @@ -17,6 +17,11 @@ namespace MapControl #endif public class Location : IEquatable { + // Arithmetic mean radius (2*a + b) / 3 == (1 - f/3) * a. + // See https://en.wikipedia.org/wiki/Earth_radius#Arithmetic_mean_radius. + // + public const double Wgs84MeanRadius = (1d - MapProjection.Wgs84Flattening / 3d) * MapProjection.Wgs84EquatorialRadius; + public Location() { } @@ -75,7 +80,7 @@ namespace MapControl /// /// Calculates great circle azimuth in degrees and distance in meters between this and the specified Location. /// - public (double, double) GetAzimuthDistance(Location location, double earthRadius = MapProjection.Wgs84MeanRadius) + public (double, double) GetAzimuthDistance(Location location, double earthRadius = Wgs84MeanRadius) { var lat1 = Latitude * Math.PI / 180d; var lon1 = Longitude * Math.PI / 180d; @@ -100,7 +105,7 @@ namespace MapControl /// /// Calculates great distance in meters between this and the specified Location. /// - public double GetDistance(Location location, double earthRadius = MapProjection.Wgs84MeanRadius) + public double GetDistance(Location location, double earthRadius = Wgs84MeanRadius) { (var _, var distance) = GetAzimuthDistance(location, earthRadius); @@ -110,7 +115,7 @@ namespace MapControl /// /// Calculates the Location on a great circle at the specified azimuth in degrees and distance in meters from this Location. /// - public Location GetLocation(double azimuth, double distance, double earthRadius = MapProjection.Wgs84MeanRadius) + public Location GetLocation(double azimuth, double distance, double earthRadius = Wgs84MeanRadius) { var lat1 = Latitude * Math.PI / 180d; var lon1 = Longitude * Math.PI / 180d; diff --git a/MapControl/Shared/MapProjection.cs b/MapControl/Shared/MapProjection.cs index 189901d8..9ee4a17e 100644 --- a/MapControl/Shared/MapProjection.cs +++ b/MapControl/Shared/MapProjection.cs @@ -31,11 +31,6 @@ namespace MapControl public const double Wgs84Flattening = 1d / 298.257223563; public const double Wgs84MeterPerDegree = Wgs84EquatorialRadius * Math.PI / 180d; - // Arithmetic mean radius (2*a + b) / 3 == (1 - f/3) * a. - // See https://en.wikipedia.org/wiki/Earth_radius#Arithmetic_mean_radius. - // - public const double Wgs84MeanRadius = (1d - Wgs84Flattening / 3d) * Wgs84EquatorialRadius; - public static MapProjectionFactory Factory { get => field ??= new MapProjectionFactory(); @@ -60,6 +55,13 @@ namespace MapControl /// public string CrsId { get; protected set; } = ""; + /// + /// The earth ellipsoid semi-major axis, or spherical earth radius, in meters. + /// + public double EquatorialRadius { get; set; } = Wgs84EquatorialRadius; + + public double MeterPerDegree => EquatorialRadius * Math.PI / 180d; + private Location center; private bool updateCenter = hasCenter; diff --git a/MapControl/Shared/OrthographicProjection.cs b/MapControl/Shared/OrthographicProjection.cs index 5a76415f..a09a7059 100644 --- a/MapControl/Shared/OrthographicProjection.cs +++ b/MapControl/Shared/OrthographicProjection.cs @@ -42,13 +42,13 @@ namespace MapControl return null; } - return new Point(EarthRadius * p.X, EarthRadius * p.Y); // p.149 (20-3/4) + return new Point(EquatorialRadius * p.X, EquatorialRadius * p.Y); // p.149 (20-3/4) } public override Location MapToLocation(double x, double y) { var rho = Math.Sqrt(x * x + y * y); - var sinC = rho / EarthRadius; // p.150 (20-19) + var sinC = rho / EquatorialRadius; // p.150 (20-19) return GetLocation(x, y, rho, sinC); } diff --git a/MapControl/Shared/PolarStereographicProjection.cs b/MapControl/Shared/PolarStereographicProjection.cs index dd19baab..ac80be70 100644 --- a/MapControl/Shared/PolarStereographicProjection.cs +++ b/MapControl/Shared/PolarStereographicProjection.cs @@ -9,8 +9,8 @@ using Avalonia; namespace MapControl { /// - /// Elliptical Polar Stereographic Projection with a given scale factor at the pole and - /// optional false easting and northing, as used by the UPS North and UPS South projections. + /// Elliptical Polar Stereographic Projection with scale factor at the pole and + /// false easting and northing, as used by the UPS North and UPS South projections. /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395), p.154-163. /// public class PolarStereographicProjection : MapProjection @@ -20,7 +20,6 @@ namespace MapControl Type = MapProjectionType.Azimuthal; } - public double EquatorialRadius { get; set; } = Wgs84EquatorialRadius; public double Flattening { get; set; } = Wgs84Flattening; public double ScaleFactor { get; set; } = 0.994; public double FalseEasting { get; set; } = 2e6; diff --git a/MapControl/Shared/StereographicProjection.cs b/MapControl/Shared/StereographicProjection.cs index 66b23f6d..ce7cab88 100644 --- a/MapControl/Shared/StereographicProjection.cs +++ b/MapControl/Shared/StereographicProjection.cs @@ -39,13 +39,13 @@ namespace MapControl var p = GetProjectedPoint(latitude, longitude); var k = 2d / (1d + p.CosC); // p.157 (21-4), k0 == 1 - return new Point(EarthRadius * k * p.X, EarthRadius * k * p.Y); // p.157 (21-2/3) + return new Point(EquatorialRadius * k * p.X, EquatorialRadius * k * p.Y); // p.157 (21-2/3) } public override Location MapToLocation(double x, double y) { var rho = Math.Sqrt(x * x + y * y); - var c = 2d * Math.Atan(rho / (2d * EarthRadius)); // p.159 (21-15), k0 == 1 + var c = 2d * Math.Atan(rho / (2d * EquatorialRadius)); // p.159 (21-15), k0 == 1 return GetLocation(x, y, rho, Math.Sin(c)); } diff --git a/MapControl/Shared/TransverseMercatorProjection.cs b/MapControl/Shared/TransverseMercatorProjection.cs index b205009b..0e341ddb 100644 --- a/MapControl/Shared/TransverseMercatorProjection.cs +++ b/MapControl/Shared/TransverseMercatorProjection.cs @@ -49,7 +49,6 @@ namespace MapControl } } - public double EquatorialRadius { get; set; } = Wgs84EquatorialRadius; public double ScaleFactor { get; set; } = 0.9996; public double CentralMeridian { get; set; } public double FalseEasting { get; set; } diff --git a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs index 881e04be..18e6ed1f 100644 --- a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs +++ b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs @@ -21,31 +21,12 @@ namespace MapControl Type = MapProjectionType.TransverseCylindrical; } + public double Flattening { get; set; } = Wgs84Flattening; public double ScaleFactor { get; set; } = 0.9996; public double CentralMeridian { get; set; } public double FalseEasting { get; set; } public double FalseNorthing { get; set; } - public double EquatorialRadius - { - get; - set - { - field = value; - M0 = MeridianDistance(LatitudeOfOrigin * Math.PI / 180d); - } - } = Wgs84EquatorialRadius; - - public double Flattening - { - get; - set - { - field = value; - M0 = MeridianDistance(LatitudeOfOrigin * Math.PI / 180d); - } - } = Wgs84Flattening; - public double LatitudeOfOrigin { get; @@ -62,11 +43,11 @@ namespace MapControl var e4 = e2 * e2; var e6 = e2 * e4; - return EquatorialRadius * - ((1d - e2 / 4d - 3d / 64d * e4 - 5d / 256d * e6) * phi - - (3d / 8d * e2 + 3d / 32d * e4 + 45d / 1024d * e6) * Math.Sin(2d * phi) + - (15d / 256d * e4 + 45d / 1024d * e6) * Math.Sin(4d * phi) - - 35d / 3072d * e6 * Math.Sin(6d * phi)); // (3-21) + return EquatorialRadius * ( + (1d - e2 / 4d - e4 * 3d / 64d - e6 * 5d / 256d) * phi - + (e2 * 3d / 8d + e4 * 3d / 32d + e6 * 45d / 1024d) * Math.Sin(2d * phi) + + (e4 * 15d / 256d + e6 * 45d / 1024d) * Math.Sin(4d * phi) - + e6 * 35d / 3072d * Math.Sin(6d * phi)); // (3-21) } public override Matrix RelativeScale(double latitude, double longitude) @@ -146,12 +127,12 @@ namespace MapControl var e14 = e1 * e13; var M = M0 + (y - FalseNorthing) / ScaleFactor; // (8-20) - var mu = M / (EquatorialRadius * (1d - e2 / 4d - 3d * e4 / 64d - 5d * e6 / 256d)); // (7-19) + var mu = M / (EquatorialRadius * (1d - e2 / 4d - e4 * 3d / 64d - e6 * 5d / 256d)); // (7-19) var phi1 = mu + - (3d * e1 / 2d - 27d * e13 / 32d) * Math.Sin(2d * mu) + - (21d * e12 / 16d - 55d * e14 / 32d) * Math.Sin(4d * mu) + - 151d * e13 / 96d * Math.Sin(6d * mu) + - 1097d * e14 / 512d * Math.Sin(8d * mu); // (3-26) + (e1 * 3d / 2d - e13 * 27d / 32d) * Math.Sin(2d * mu) + + (e12 * 21d / 16d - e14 * 55d / 32d) * Math.Sin(4d * mu) + + e13 * 151d / 96d * Math.Sin(6d * mu) + + e14 * 1097d / 512d * Math.Sin(8d * mu); // (3-26) var sinPhi1 = Math.Sin(phi1); var cosPhi1 = Math.Cos(phi1); diff --git a/MapControl/Shared/WebMercatorProjection.cs b/MapControl/Shared/WebMercatorProjection.cs index 7faa3dd5..f86cfe4e 100644 --- a/MapControl/Shared/WebMercatorProjection.cs +++ b/MapControl/Shared/WebMercatorProjection.cs @@ -37,15 +37,15 @@ namespace MapControl public override Point? LocationToMap(double latitude, double longitude) { return new Point( - Wgs84MeterPerDegree * longitude, - Wgs84MeterPerDegree * LatitudeToY(latitude)); + MeterPerDegree * longitude, + MeterPerDegree * LatitudeToY(latitude)); } public override Location MapToLocation(double x, double y) { return new Location( - YToLatitude(y / Wgs84MeterPerDegree), - x / Wgs84MeterPerDegree); + YToLatitude(y / MeterPerDegree), + x / MeterPerDegree); } public static double LatitudeToY(double latitude) diff --git a/MapControl/Shared/WorldMercatorProjection.cs b/MapControl/Shared/WorldMercatorProjection.cs index b4bfee85..77b645c8 100644 --- a/MapControl/Shared/WorldMercatorProjection.cs +++ b/MapControl/Shared/WorldMercatorProjection.cs @@ -39,15 +39,15 @@ namespace MapControl public override Point? LocationToMap(double latitude, double longitude) { return new Point( - Wgs84MeterPerDegree * longitude, - Wgs84MeterPerDegree * LatitudeToY(latitude)); + MeterPerDegree * longitude, + MeterPerDegree * LatitudeToY(latitude)); } public override Location MapToLocation(double x, double y) { return new Location( - YToLatitude(y / Wgs84MeterPerDegree), - x / Wgs84MeterPerDegree); + YToLatitude(y / MeterPerDegree), + x / MeterPerDegree); } public static double RelativeScale(double latitude) @@ -93,10 +93,10 @@ namespace MapControl var chi = Math.PI / 2d - 2d * Math.Atan(t); // p.45 (7-13) return chi + - (e2 / 2d + 5d * e4 / 24d + e6 / 12d + 13d * e8 / 360d) * Math.Sin(2d * chi) + - (7d * e4 / 48d + 29d * e6 / 240d + 811d * e8 / 11520d) * Math.Sin(4d * chi) + - (7d * e6 / 120d + 81d * e8 / 1120d) * Math.Sin(6d * chi) + - 4279d * e8 / 161280d * Math.Sin(8d * chi); // p.45 (3-5) + (e2 / 2d + e4 * 5d / 24d + e6 / 12d + e8 * 13d / 360d) * Math.Sin(2d * chi) + + (e4 * 7d / 48d + e6 * 29d / 240d + e8 * 811d / 11520d) * Math.Sin(4d * chi) + + (e6 * 7d / 120d + e8 * 81d / 1120d) * Math.Sin(6d * chi) + + e8 * 4279d / 161280d * Math.Sin(8d * chi); // p.45 (3-5) } } } diff --git a/SampleApps/WpfApplication/MainWindow.xaml b/SampleApps/WpfApplication/MainWindow.xaml index 6ab76a16..4293f0a5 100644 --- a/SampleApps/WpfApplication/MainWindow.xaml +++ b/SampleApps/WpfApplication/MainWindow.xaml @@ -211,8 +211,10 @@ Description="© [BKG](https://gdz.bkg.bund.de/index.php/default/webdienste/topplus-produkte/wms-topplusopen-mit-layer-fur-normalausgabe-und-druck-wms-topplus-open.html)"/> +