diff --git a/MapControl/Shared/AzimuthalEquidistantProjection.cs b/MapControl/Shared/AzimuthalEquidistantProjection.cs index a403d396..8f1fc91a 100644 --- a/MapControl/Shared/AzimuthalEquidistantProjection.cs +++ b/MapControl/Shared/AzimuthalEquidistantProjection.cs @@ -11,7 +11,7 @@ namespace MapControl /// Spherical Azimuthal Equidistant Projection - No standard CRS identifier. /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395), p.195-197. /// - public class AzimuthalEquidistantProjection : MapProjection + public class AzimuthalEquidistantProjection : AzimuthalProjection { public const string DefaultCrsId = "AUTO2:97003"; // proprietary CRS identifier @@ -22,38 +22,43 @@ namespace MapControl public AzimuthalEquidistantProjection(string crsId) { - Type = MapProjectionType.Azimuthal; CrsId = crsId; } - public double EarthRadius { get; set; } = Wgs84MeanRadius; + public override Point RelativeScale(double latitude, double longitude) + { + (var cosC, var _, var _) = GetPointValues(latitude, longitude); + var k = 1d; + + if (cosC < 1d) + { + var c = Math.Acos(cosC); + k = c / Math.Sin(c); // p.195 (25-2) + } + + return new Point(k, k); + } public override Point? LocationToMap(double latitude, double longitude) { - if (Center.Equals(latitude, longitude)) + (var cosC, var x, var y) = GetPointValues(latitude, longitude); + var k = 1d; + + if (cosC < 1d) { - return new Point(); + var c = Math.Acos(cosC); + k = c / Math.Sin(c); // p.195 (25-2) } - Center.GetAzimuthDistance(latitude, longitude, out double azimuth, out double distance); - - var mapDistance = distance * EarthRadius; - - return new Point(mapDistance * Math.Sin(azimuth), mapDistance * Math.Cos(azimuth)); + return new Point(EarthRadius * k * x, EarthRadius * k * y); // p.195 (22-4/5) } public override Location MapToLocation(double x, double y) { - if (x == 0d && y == 0d) - { - return new Location(Center.Latitude, Center.Longitude); - } + var rho = Math.Sqrt(x * x + y * y); + var c = rho / EarthRadius; // p.196 (25-15) - var azimuth = Math.Atan2(x, y); - var mapDistance = Math.Sqrt(x * x + y * y); - var distance = mapDistance / EarthRadius; - - return Center.GetLocation(azimuth, distance); + return GetLocation(x, y, rho, Math.Sin(c)); } } } diff --git a/MapControl/Shared/AzimuthalProjection.cs b/MapControl/Shared/AzimuthalProjection.cs new file mode 100644 index 00000000..825f88fd --- /dev/null +++ b/MapControl/Shared/AzimuthalProjection.cs @@ -0,0 +1,47 @@ +using System; + +namespace MapControl +{ + /// + /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395). + /// + public abstract class AzimuthalProjection : MapProjection + { + protected AzimuthalProjection() + { + Type = MapProjectionType.Azimuthal; + } + + public double EarthRadius { get; set; } = Wgs84MeanRadius; + + protected (double, double, double) GetPointValues(double latitude, double longitude) + { + 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; + + return (cosC, x, y); + } + + protected Location GetLocation(double x, double y, double rho, double sinC) + { + var cosC = Math.Sqrt(1d - sinC * sinC); + var phi1 = Center.Latitude * Math.PI / 180d; + var cosPhi1 = Math.Cos(phi1); + var sinPhi1 = Math.Sin(phi1); + var phi = Math.Asin(cosC * sinPhi1 + y * sinC * cosPhi1 / rho); // (20-14) + var dLambda = Math.Atan2(x * sinC, rho * cosPhi1 * cosC - y * sinPhi1 * sinC); // (20-15), λ - λ0 + + return new Location(phi * 180d / Math.PI, dLambda * 180d / Math.PI + Center.Longitude); + } + } +} diff --git a/MapControl/Shared/GnomonicProjection.cs b/MapControl/Shared/GnomonicProjection.cs index 3f21fa51..46399eb7 100644 --- a/MapControl/Shared/GnomonicProjection.cs +++ b/MapControl/Shared/GnomonicProjection.cs @@ -11,7 +11,7 @@ namespace MapControl /// Spherical Gnomonic Projection - AUTO2:97001. /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395), p.165-167. /// - public class GnomonicProjection : MapProjection + public class GnomonicProjection : AzimuthalProjection { public const string DefaultCrsId = "AUTO2:97001"; // GeoServer non-standard CRS identifier @@ -22,40 +22,37 @@ namespace MapControl public GnomonicProjection(string crsId) { - Type = MapProjectionType.Azimuthal; CrsId = crsId; } - public double EarthRadius { get; set; } = Wgs84MeanRadius; + public override Point RelativeScale(double latitude, double longitude) + { + (var cosC, var _, var _) = GetPointValues(latitude, longitude); + var h = 1d / (cosC * cosC); // p.165 (22-2) + + return new Point(h, h); // TODO: rotate + } public override Point? LocationToMap(double latitude, double longitude) { - if (Center.Equals(latitude, longitude)) + (var cosC, var x, var y) = GetPointValues(latitude, longitude); + + if (cosC <= 0d) { - return new Point(); + return null; } - Center.GetAzimuthDistance(latitude, longitude, out double azimuth, out double distance); + var k = 1d / cosC; // p.165 (22-3) - var mapDistance = distance < Math.PI / 2d - ? Math.Tan(distance) * EarthRadius - : double.PositiveInfinity; - - return new Point(mapDistance * Math.Sin(azimuth), mapDistance * Math.Cos(azimuth)); + return new Point(EarthRadius * k * x, EarthRadius * k * y); // p.165 (22-4/5) } public override Location MapToLocation(double x, double y) { - if (x == 0d && y == 0d) - { - return new Location(Center.Latitude, Center.Longitude); - } + var rho = Math.Sqrt(x * x + y * y); + var c = Math.Atan(rho / EarthRadius); // p.167 (22-16) - var azimuth = Math.Atan2(x, y); - var mapDistance = Math.Sqrt(x * x + y * y); - var distance = Math.Atan(mapDistance / EarthRadius); - - return Center.GetLocation(azimuth, distance); + return GetLocation(x, y, rho, Math.Sin(c)); } } } diff --git a/MapControl/Shared/OrthographicProjection.cs b/MapControl/Shared/OrthographicProjection.cs index cbd7e487..10181779 100644 --- a/MapControl/Shared/OrthographicProjection.cs +++ b/MapControl/Shared/OrthographicProjection.cs @@ -11,7 +11,7 @@ namespace MapControl /// Spherical Orthographic Projection - AUTO2:42003. /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395), p.148-150. /// - public class OrthographicProjection : MapProjection + public class OrthographicProjection : AzimuthalProjection { public const string DefaultCrsId = "AUTO2:42003"; @@ -22,62 +22,30 @@ namespace MapControl public OrthographicProjection(string crsId) { - Type = MapProjectionType.Azimuthal; CrsId = crsId; } - public double EarthRadius { get; set; } = Wgs84MeanRadius; + public override Point RelativeScale(double latitude, double longitude) + { + (var cosC, var x, var y) = GetPointValues(latitude, longitude); + var h = cosC; // p.149 (20-5) + + return new Point(h, h); + } public override Point? LocationToMap(double latitude, double longitude) { - if (Center.Equals(latitude, longitude)) - { - return new Point(); - } + (var cosC, var x, var y) = GetPointValues(latitude, longitude); - var phi = latitude * Math.PI / 180d; - var phi1 = Center.Latitude * Math.PI / 180d; - var dLambda = (longitude - Center.Longitude) * Math.PI / 180d; // λ - λ0 - - if (Math.Abs(phi - phi1) > Math.PI / 2d || Math.Abs(dLambda) > Math.PI / 2d) - { - return null; - } - - var x = EarthRadius * Math.Cos(phi) * Math.Sin(dLambda); // p.149 (20-3) - var y = EarthRadius * (Math.Cos(phi1) * Math.Sin(phi) - - Math.Sin(phi1) * Math.Cos(phi) * Math.Cos(dLambda)); // p.149 (20-4) - return new Point(x, y); + return cosC >= 0d ? new Point(EarthRadius * x, EarthRadius * y) : null; // p.149 (20-3/4) } public override Location MapToLocation(double x, double y) { - if (x == 0d && y == 0d) - { - return new Location(Center.Latitude, Center.Longitude); - } + var rho = Math.Sqrt(x * x + y * y); + var sinC = rho / EarthRadius; // p.150 (20-19) - x /= EarthRadius; - y /= EarthRadius; - var r2 = x * x + y * y; - - if (r2 > 1d) - { - return null; - } - - var r = Math.Sqrt(r2); // p.150 (20-18), r=ρ/R - var sinC = r; // p.150 (20-19) - var cosC = Math.Sqrt(1d - r2); - - var phi1 = Center.Latitude * Math.PI / 180d; - var cosPhi1 = Math.Cos(phi1); - var sinPhi1 = Math.Sin(phi1); - - var phi = Math.Asin(cosC * sinPhi1 + y * sinC * cosPhi1 / r); // p.150 (20-14) - var dLambda = Math.Atan2(x * sinC, r * cosC * cosPhi1 - y * sinC * sinPhi1); // p.150 (20-15) - - return new Location(180d / Math.PI * phi, 180d / Math.PI * dLambda + Center.Longitude); + return sinC <= 1d ? GetLocation(x, y, rho, sinC) : null; } } } diff --git a/MapControl/Shared/StereographicProjection.cs b/MapControl/Shared/StereographicProjection.cs index b573398a..473e20f6 100644 --- a/MapControl/Shared/StereographicProjection.cs +++ b/MapControl/Shared/StereographicProjection.cs @@ -11,7 +11,7 @@ namespace MapControl /// Spherical Stereographic Projection - AUTO2:97002. /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395), p.157-160. /// - public class StereographicProjection : MapProjection + public class StereographicProjection : AzimuthalProjection { public const string DefaultCrsId = "AUTO2:97002"; // GeoServer non-standard CRS identifier @@ -22,38 +22,31 @@ namespace MapControl public StereographicProjection(string crsId) { - Type = MapProjectionType.Azimuthal; CrsId = crsId; } - public double EarthRadius { get; set; } = Wgs84MeanRadius; + public override Point RelativeScale(double latitude, double longitude) + { + (var cosC, var _, var _) = GetPointValues(latitude, longitude); + var k = 2d / (1d + cosC); // p.157 (21-4), k0 == 1 + + return new Point(k, k); + } public override Point? LocationToMap(double latitude, double longitude) { - if (Center.Equals(latitude, longitude)) - { - return new Point(); - } + (var cosC, var x, var y) = GetPointValues(latitude, longitude); + var k = 2d / (1d + cosC); // p.157 (21-4), k0 == 1 - Center.GetAzimuthDistance(latitude, longitude, out double azimuth, out double distance); - - var mapDistance = Math.Tan(distance / 2d) * 2d * EarthRadius; - - return new Point(mapDistance * Math.Sin(azimuth), mapDistance * Math.Cos(azimuth)); + return new Point(EarthRadius * k * x, EarthRadius * k * y); // p.157 (21-2/3) } public override Location MapToLocation(double x, double y) { - if (x == 0d && y == 0d) - { - return new Location(Center.Latitude, Center.Longitude); - } + var rho = Math.Sqrt(x * x + y * y); + var c = 2d * Math.Atan(rho / (2d * EarthRadius)); // p.159 (21-15), k0 == 1 - var azimuth = Math.Atan2(x, y); - var mapDistance = Math.Sqrt(x * x + y * y); - var distance = 2d * Math.Atan(mapDistance / (2d * EarthRadius)); - - return Center.GetLocation(azimuth, distance); + return GetLocation(x, y, rho, Math.Sin(c)); } } } diff --git a/MapControl/Shared/TransverseMercatorProjection.cs b/MapControl/Shared/TransverseMercatorProjection.cs index 5c5488d3..2980de5e 100644 --- a/MapControl/Shared/TransverseMercatorProjection.cs +++ b/MapControl/Shared/TransverseMercatorProjection.cs @@ -83,7 +83,7 @@ namespace MapControl // λ - λ0 var dLambda = (longitude - CentralMeridian) * Math.PI / 180d; // ξ' - var xi_ = Math.Atan(t / Math.Cos(dLambda)); + var xi_ = Math.Atan2(t, Math.Cos(dLambda)); // η' var eta_ = Atanh(Math.Sin(dLambda) / Math.Sqrt(1d + t * t)); // ξ @@ -128,7 +128,7 @@ namespace MapControl d2 * Math.Sin(4d * chi) + d3 * Math.Sin(6d * chi); // λ - λ0 - var dLambda = Math.Atan(Math.Sinh(eta_) / Math.Cos(xi_)); + var dLambda = Math.Atan2(Math.Sinh(eta_), Math.Cos(xi_)); return new Location( phi * 180d / Math.PI,