diff --git a/MapControl/Shared/PolarStereographicProjection.cs b/MapControl/Shared/PolarStereographicProjection.cs index 1c436ea7..ff1b4851 100644 --- a/MapControl/Shared/PolarStereographicProjection.cs +++ b/MapControl/Shared/PolarStereographicProjection.cs @@ -69,12 +69,13 @@ namespace MapControl x = sign * (x - FalseEasting); y = sign * (y - FalseNorthing); - var e = Math.Sqrt((2d - Flattening) * Flattening); + var e2 = (2d - Flattening) * Flattening; + var e = Math.Sqrt(e2); var r = Math.Sqrt(x * x + y * y); // p.162 (20-18) var t = r * Math.Sqrt(Math.Pow(1d + e, 1d + e) * Math.Pow(1d - e, 1d - e)) / (2d * EquatorialRadius * ScaleFactor); // p.162 (21-39) - var phi = WorldMercatorProjection.ApproximateLatitude(e, t); // p.162 (3-5) + var phi = WorldMercatorProjection.ApproximateLatitude(e2, t); // p.162 (3-5) var lambda = Math.Atan2(x, -y); // p.162 (20-16) return new Location(sign * phi * 180d / Math.PI, sign * lambda * 180d / Math.PI); diff --git a/MapControl/Shared/WebMercatorProjection.cs b/MapControl/Shared/WebMercatorProjection.cs index 12f68f4b..f1e1db50 100644 --- a/MapControl/Shared/WebMercatorProjection.cs +++ b/MapControl/Shared/WebMercatorProjection.cs @@ -43,8 +43,8 @@ namespace MapControl public override Location MapToLocation(double x, double y) { - return new Location( - YToLatitude(y / EquatorialRadius * 180d / Math.PI), + return new Location(YToLatitude( + y / EquatorialRadius * 180d / Math.PI), x / EquatorialRadius * 180d / Math.PI); } diff --git a/MapControl/Shared/WorldMercatorProjection.cs b/MapControl/Shared/WorldMercatorProjection.cs index 0d03ad16..631f5bb7 100644 --- a/MapControl/Shared/WorldMercatorProjection.cs +++ b/MapControl/Shared/WorldMercatorProjection.cs @@ -14,10 +14,10 @@ namespace MapControl /// public class WorldMercatorProjection : MapProjection { - public static readonly double Wgs84Eccentricity = Math.Sqrt((2d - Wgs84Flattening) * Wgs84Flattening); - public const string DefaultCrsId = "EPSG:3395"; + public double Flattening { get; set; } = Wgs84Flattening; + public WorldMercatorProjection() // parameterless constructor for XAML : this(DefaultCrsId) { @@ -31,62 +31,51 @@ namespace MapControl public override Matrix RelativeTransform(double latitude, double longitude) { - var k = RelativeScale(latitude); + var phi = latitude * Math.PI / 180d; + var e2 = (2d - Flattening) * Flattening; + var sinPhi = Math.Sin(phi); + var k = Math.Sqrt(1d - e2 * sinPhi * sinPhi) / Math.Cos(phi); // p.44 (7-8) return new Matrix(k, 0d, 0d, k, 0d, 0d); } public override Point LocationToMap(double latitude, double longitude) { - return new Point( - EquatorialRadius * Math.PI / 180d * longitude, - EquatorialRadius * Math.PI / 180d * LatitudeToY(latitude)); + var x = EquatorialRadius * longitude * Math.PI / 180d; + double y; + + if (latitude <= -90d) + { + y = double.NegativeInfinity; + } + else if (latitude >= 90d) + { + y = double.PositiveInfinity; + } + else + { + var phi = latitude * Math.PI / 180d; + var e = Math.Sqrt((2d - Flattening) * Flattening); + var eSinPhi = e * Math.Sin(phi); + var f = Math.Pow((1d - eSinPhi) / (1d + eSinPhi), e / 2d); + + y = EquatorialRadius * Math.Log(Math.Tan(phi / 2d + Math.PI / 4d) * f); // p.44 (7-7) + } + + return new Point(x, y); } public override Location MapToLocation(double x, double y) { - return new Location( - YToLatitude(y / EquatorialRadius * 180d / Math.PI), - x / EquatorialRadius * 180d / Math.PI); + var t = Math.Exp(-y / EquatorialRadius); // p.44 (7-10) + var phi = ApproximateLatitude((2d - Flattening) * Flattening, t); // p.45 (3-5) + var lambda = x / EquatorialRadius; + + return new Location(phi * 180d / Math.PI, lambda * 180d / Math.PI); } - public static double RelativeScale(double latitude) + internal static double ApproximateLatitude(double e2, double t) { - var phi = latitude * Math.PI / 180d; - var eSinPhi = Wgs84Eccentricity * Math.Sin(phi); - - return Math.Sqrt(1d - eSinPhi * eSinPhi) / Math.Cos(phi); // p.44 (7-8) - } - - public static double LatitudeToY(double latitude) - { - if (latitude <= -90d) - { - return double.NegativeInfinity; - } - - if (latitude >= 90d) - { - return double.PositiveInfinity; - } - - var phi = latitude * Math.PI / 180d; - var eSinPhi = Wgs84Eccentricity * Math.Sin(phi); - var f = Math.Pow((1d - eSinPhi) / (1d + eSinPhi), Wgs84Eccentricity / 2d); - - return Math.Log(Math.Tan(phi / 2d + Math.PI / 4d) * f) * 180d / Math.PI; // p.44 (7-7) - } - - public static double YToLatitude(double y) - { - var t = Math.Exp(-y * Math.PI / 180d); // p.44 (7-10) - - return ApproximateLatitude(Wgs84Eccentricity, t) * 180d / Math.PI; - } - - internal static double ApproximateLatitude(double e, double t) - { - var e2 = e * e; var e4 = e2 * e2; var e6 = e2 * e4; var e8 = e2 * e6; diff --git a/MapProjections/Shared/ProjNetMapProjection.cs b/MapProjections/Shared/ProjNetMapProjection.cs index 9c813939..946aeba1 100644 --- a/MapProjections/Shared/ProjNetMapProjection.cs +++ b/MapProjections/Shared/ProjNetMapProjection.cs @@ -15,7 +15,7 @@ namespace MapControl.Projections /// public class ProjNetMapProjection : MapProjection { - protected MapProjection FallbackProjection { get; } + protected MapProjection FallbackProjection { get; private set; } protected ProjNetMapProjection(MapProjection fallbackProjection) { @@ -49,10 +49,10 @@ namespace MapControl.Projections field = value ?? throw new ArgumentNullException(nameof(value)); - var name = field.Projection?.Name ?? + var projection = field.Projection ?? throw new ArgumentException("CoordinateSystem.Projection must not be null.", nameof(value)); - IsNormalCylindrical = name.StartsWith("Mercator") || name.Contains("Pseudo-Mercator"); + IsNormalCylindrical = projection.Name.StartsWith("Mercator") || projection.Name.Contains("Pseudo-Mercator"); var transformFactory = new CoordinateTransformationFactory(); @@ -67,6 +67,48 @@ namespace MapControl.Projections CrsId = !string.IsNullOrEmpty(field.Authority) && field.AuthorityCode > 0 ? $"{field.Authority}:{field.AuthorityCode}" : string.Empty; + + var ellipsoid = field.GeographicCoordinateSystem.HorizontalDatum.Ellipsoid; + + if (projection.Name.Contains("Pseudo-Mercator")) + { + FallbackProjection = new MapControl.WebMercatorProjection + { + EquatorialRadius = ellipsoid.SemiMajorAxis + }; + } + else if (projection.Name.StartsWith("Mercator")) + { + FallbackProjection = new MapControl.WorldMercatorProjection + { + EquatorialRadius = ellipsoid.SemiMajorAxis, + Flattening = 1d / ellipsoid.InverseFlattening + }; + } + else if (projection.Name.StartsWith("Transverse_Mercator")) + { + FallbackProjection = new TransverseMercatorProjection + { + EquatorialRadius = ellipsoid.SemiMajorAxis, + Flattening = 1d / ellipsoid.InverseFlattening, + CentralMeridian = projection.GetParameter("central_meridian").Value, + ScaleFactor = projection.GetParameter("scale_factor").Value, + FalseEasting = projection.GetParameter("false_easting").Value, + FalseNorthing = projection.GetParameter("false_northing").Value + }; + } + else if (projection.Name.StartsWith("Polar_Stereographic")) + { + FallbackProjection = new PolarStereographicProjection + { + EquatorialRadius = ellipsoid.SemiMajorAxis, + Flattening = 1d / ellipsoid.InverseFlattening, + ScaleFactor = projection.GetParameter("scale_factor").Value, + FalseEasting = projection.GetParameter("false_easting").Value, + FalseNorthing = projection.GetParameter("false_northing").Value, + Hemisphere = projection.GetParameter("latitude_of_origin").Value >= 0 ? Hemisphere.North : Hemisphere.South + }; + } } }