Updated map projections

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);

MapControl/Shared/WebMercatorProjection.cs

@@ -43,8 +43,8 @@ namespace MapControl
 
         public override Location MapToLocation(double x, double y)
         {
-            return new Location(
+            return new Location(YToLatitude(
-                YToLatitude(y / EquatorialRadius * 180d / Math.PI),
+                y / EquatorialRadius * 180d / Math.PI),
                 x / EquatorialRadius * 180d / Math.PI);
         }
 

MapControl/Shared/WorldMercatorProjection.cs

@@ -14,10 +14,10 @@ namespace MapControl
     /// </summary>
     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(
+            var x = EquatorialRadius * longitude * Math.PI / 180d;
-                EquatorialRadius * Math.PI / 180d * longitude,
+            double y;
-                EquatorialRadius * Math.PI / 180d * LatitudeToY(latitude));
+
+            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(
+            var t = Math.Exp(-y / EquatorialRadius); // p.44 (7-10)
-                YToLatitude(y / EquatorialRadius * 180d / Math.PI),
+            var phi = ApproximateLatitude((2d - Flattening) * Flattening, t); // p.45 (3-5)
-                x / EquatorialRadius * 180d / Math.PI);
+            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;

MapProjections/Shared/ProjNetMapProjection.cs

@@ -15,7 +15,7 @@ namespace MapControl.Projections
     /// </summary>
     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
+                    };
+                }
             }
         }