diff --git a/MapControl/Shared/AzimuthalEquidistantProjection.cs b/MapControl/Shared/AzimuthalEquidistantProjection.cs index 8e522514..1e6db43f 100644 --- a/MapControl/Shared/AzimuthalEquidistantProjection.cs +++ b/MapControl/Shared/AzimuthalEquidistantProjection.cs @@ -28,12 +28,12 @@ namespace MapControl public override Matrix RelativeScale(double latitude, double longitude) { - (var cosC, var _, var _) = GetPointValues(latitude, longitude); + var p = GetProjectedPoint(latitude, longitude); var k = 1d; - if (cosC < 1d) + if (p.CosC < 1d) { - var c = Math.Acos(cosC); + var c = Math.Acos(p.CosC); k = c / Math.Sin(c); // p.195 (25-2) } @@ -42,16 +42,16 @@ namespace MapControl public override Point? LocationToMap(double latitude, double longitude) { - (var cosC, var x, var y) = GetPointValues(latitude, longitude); + var p = GetProjectedPoint(latitude, longitude); var k = 1d; - if (cosC < 1d) + if (p.CosC < 1d) { - var c = Math.Acos(cosC); + var c = Math.Acos(p.CosC); k = c / Math.Sin(c); // p.195 (25-2) } - return new Point(EarthRadius * k * x, EarthRadius * k * y); // p.195 (22-4/5) + return new Point(EarthRadius * k * p.X, EarthRadius * k * p.Y); // p.195 (22-4/5) } public override Location MapToLocation(double x, double y) diff --git a/MapControl/Shared/AzimuthalProjection.cs b/MapControl/Shared/AzimuthalProjection.cs index 825f88fd..c24bd069 100644 --- a/MapControl/Shared/AzimuthalProjection.cs +++ b/MapControl/Shared/AzimuthalProjection.cs @@ -3,7 +3,7 @@ namespace MapControl { /// - /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395). + /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395), p.141. /// public abstract class AzimuthalProjection : MapProjection { @@ -14,22 +14,33 @@ namespace MapControl public double EarthRadius { get; set; } = Wgs84MeanRadius; - protected (double, double, double) GetPointValues(double latitude, double longitude) + public readonly struct ProjectedPoint { - var phi = latitude * Math.PI / 180d; - var phi1 = Center.Latitude * Math.PI / 180d; - var dLambda = (longitude - Center.Longitude) * Math.PI / 180d; // λ - λ0 - var cosPhi = Math.Cos(phi); - var sinPhi = Math.Sin(phi); - var cosPhi1 = Math.Cos(phi1); - var sinPhi1 = Math.Sin(phi1); - var cosLambda = Math.Cos(dLambda); - var sinLambda = Math.Sin(dLambda); - var cosC = sinPhi1 * sinPhi + cosPhi1 * cosPhi * cosLambda; // (5-3) - var x = cosPhi * sinLambda; - var y = cosPhi1 * sinPhi - sinPhi1 * cosPhi * cosLambda; + public double X { get; } + public double Y { get; } + public double CosC { get; } - return (cosC, x, y); + public ProjectedPoint(double centerLatitude, double centerLongitude, double latitude, double longitude) + { + var phi = latitude * Math.PI / 180d; + var phi1 = centerLatitude * Math.PI / 180d; + var dLambda = (longitude - centerLongitude) * Math.PI / 180d; // λ - λ0 + var cosPhi = Math.Cos(phi); + var sinPhi = Math.Sin(phi); + var cosPhi1 = Math.Cos(phi1); + var sinPhi1 = Math.Sin(phi1); + var cosLambda = Math.Cos(dLambda); + var sinLambda = Math.Sin(dLambda); + + X = cosPhi * sinLambda; + Y = cosPhi1 * sinPhi - sinPhi1 * cosPhi * cosLambda; + CosC = sinPhi1 * sinPhi + cosPhi1 * cosPhi * cosLambda; // (5-3) + } + } + + protected ProjectedPoint GetProjectedPoint(double latitude, double longitude) + { + return new ProjectedPoint(Center.Latitude, Center.Longitude, latitude, longitude); } protected Location GetLocation(double x, double y, double rho, double sinC) diff --git a/MapControl/Shared/GnomonicProjection.cs b/MapControl/Shared/GnomonicProjection.cs index b4c83a6d..0f264b3e 100644 --- a/MapControl/Shared/GnomonicProjection.cs +++ b/MapControl/Shared/GnomonicProjection.cs @@ -28,27 +28,27 @@ namespace MapControl public override Matrix RelativeScale(double latitude, double longitude) { - (var cosC, var x, var y) = GetPointValues(latitude, longitude); - var k = 1d / cosC; // p.165 (22-3) + var p = GetProjectedPoint(latitude, longitude); + var k = 1d / p.CosC; // p.165 (22-3) var h = k * k; // p.165 (22-2) var scale = new Matrix(h, 0d, 0d, k, 0d, 0d); - scale.Rotate(-Math.Atan2(y, x) * 180d / Math.PI); + scale.Rotate(-Math.Atan2(p.Y, p.X) * 180d / Math.PI); return scale; } public override Point? LocationToMap(double latitude, double longitude) { - (var cosC, var x, var y) = GetPointValues(latitude, longitude); + var p = GetProjectedPoint(latitude, longitude); - if (cosC <= 0d) + if (p.CosC <= 0d) { return null; } - var k = 1d / cosC; // p.165 (22-3) + var k = 1d / p.CosC; // p.165 (22-3) - return new Point(EarthRadius * k * x, EarthRadius * k * y); // p.165 (22-4/5) + return new Point(EarthRadius * k * p.X, EarthRadius * k * p.Y); // p.165 (22-4/5) } public override Location MapToLocation(double x, double y) diff --git a/MapControl/Shared/MapProjection.cs b/MapControl/Shared/MapProjection.cs index 8657486f..3e9f53eb 100644 --- a/MapControl/Shared/MapProjection.cs +++ b/MapControl/Shared/MapProjection.cs @@ -28,8 +28,8 @@ namespace MapControl public abstract class MapProjection { public const double Wgs84EquatorialRadius = 6378137d; - public const double Wgs84MeterPerDegree = Wgs84EquatorialRadius * Math.PI / 180d; 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. diff --git a/MapControl/Shared/OrthographicProjection.cs b/MapControl/Shared/OrthographicProjection.cs index c94d425d..3c7e401c 100644 --- a/MapControl/Shared/OrthographicProjection.cs +++ b/MapControl/Shared/OrthographicProjection.cs @@ -28,19 +28,19 @@ namespace MapControl public override Matrix RelativeScale(double latitude, double longitude) { - (var cosC, var x, var y) = GetPointValues(latitude, longitude); - var h = cosC; // p.149 (20-5) + var p = GetProjectedPoint(latitude, longitude); + var h = p.CosC; // p.149 (20-5) var scale = new Matrix(h, 0d, 0d, 1d, 0d, 0d); - scale.Rotate(-Math.Atan2(y, x) * 180d / Math.PI); + scale.Rotate(-Math.Atan2(p.Y, p.X) * 180d / Math.PI); return scale; } public override Point? LocationToMap(double latitude, double longitude) { - (var cosC, var x, var y) = GetPointValues(latitude, longitude); + var p = GetProjectedPoint(latitude, longitude); - return cosC >= 0d ? new Point(EarthRadius * x, EarthRadius * y) : null; // p.149 (20-3/4) + return p.CosC >= 0d ? new Point(EarthRadius * p.X, EarthRadius * p.Y) : null; // p.149 (20-3/4) } public override Location MapToLocation(double x, double y) diff --git a/MapControl/Shared/StereographicProjection.cs b/MapControl/Shared/StereographicProjection.cs index a1304be7..b166c489 100644 --- a/MapControl/Shared/StereographicProjection.cs +++ b/MapControl/Shared/StereographicProjection.cs @@ -28,18 +28,18 @@ namespace MapControl public override Matrix RelativeScale(double latitude, double longitude) { - (var cosC, var _, var _) = GetPointValues(latitude, longitude); - var k = 2d / (1d + cosC); // p.157 (21-4), k0 == 1 + var p = GetProjectedPoint(latitude, longitude); + var k = 2d / (1d + p.CosC); // p.157 (21-4), k0 == 1 return new Matrix(k, 0d, 0d, k, 0d, 0d); } public override Point? LocationToMap(double latitude, double longitude) { - (var cosC, var x, var y) = GetPointValues(latitude, longitude); - var k = 2d / (1d + cosC); // p.157 (21-4), k0 == 1 + var p = GetProjectedPoint(latitude, longitude); + var k = 2d / (1d + p.CosC); // p.157 (21-4), k0 == 1 - return new Point(EarthRadius * k * x, EarthRadius * k * y); // p.157 (21-2/3) + return new Point(EarthRadius * k * p.X, EarthRadius * k * p.Y); // p.157 (21-2/3) } public override Location MapToLocation(double x, double y) diff --git a/MapProjections/Shared/Wgs84OrthographicProjection.cs b/MapProjections/Shared/Wgs84OrthographicProjection.cs index 57458f59..1a22b650 100644 --- a/MapProjections/Shared/Wgs84OrthographicProjection.cs +++ b/MapProjections/Shared/Wgs84OrthographicProjection.cs @@ -1,7 +1,15 @@ -using System.Globalization; +using System; +using System.Globalization; +#if WPF +using System.Windows.Media; +#endif namespace MapControl.Projections { + /// + /// Spherical Orthographic Projection - AUTO2:42003. + /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395), p.148-150. + /// public class Wgs84OrthographicProjection : ProjNetMapProjection { public Wgs84OrthographicProjection() @@ -25,5 +33,15 @@ namespace MapControl.Projections CoordinateSystemWkt = string.Format( CultureInfo.InvariantCulture, wktFormat, Center.Latitude, Center.Longitude); } + + public override Matrix RelativeScale(double latitude, double longitude) + { + var p = new AzimuthalProjection.ProjectedPoint(Center.Latitude, Center.Longitude, latitude, longitude); + var h = p.CosC; // p.149 (20-5) + + var scale = new Matrix(h, 0d, 0d, 1d, 0d, 0d); + scale.Rotate(-Math.Atan2(p.Y, p.X) * 180d / Math.PI); + return scale; + } } } diff --git a/MapProjections/Shared/Wgs84StereographicProjection.cs b/MapProjections/Shared/Wgs84StereographicProjection.cs index 0ed274df..54408fcf 100644 --- a/MapProjections/Shared/Wgs84StereographicProjection.cs +++ b/MapProjections/Shared/Wgs84StereographicProjection.cs @@ -1,7 +1,14 @@ using System.Globalization; +#if WPF +using System.Windows.Media; +#endif namespace MapControl.Projections { + /// + /// Spherical Stereographic Projection - AUTO2:97002. + /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395), p.157-160. + /// public class Wgs84StereographicProjection : ProjNetMapProjection { public Wgs84StereographicProjection() @@ -25,5 +32,13 @@ namespace MapControl.Projections CoordinateSystemWkt = string.Format( CultureInfo.InvariantCulture, wktFormat, Center.Latitude, Center.Longitude); } + + public override Matrix RelativeScale(double latitude, double longitude) + { + var p = new AzimuthalProjection.ProjectedPoint(Center.Latitude, Center.Longitude, latitude, longitude); + var k = 2d / (1d + p.CosC); // p.157 (21-4), k0 == 1 + + return new Matrix(k, 0d, 0d, k, 0d, 0d); + } } }