New TransverseMercatorProjection implementation

MapControl/Shared/Etrs89UtmProjection.cs

@@ -0,0 +1,40 @@
+// XAML Map Control - https://github.com/ClemensFischer/XAML-Map-Control
+// Copyright © 2024 Clemens Fischer
+// Licensed under the Microsoft Public License (Ms-PL)
+
+using System;
+
+namespace MapControl
+{
+    /// <summary>
+    /// ETRS89 UTM Projection with zone number.
+    /// </summary>
+    public class Etrs89UtmProjection : TransverseMercatorProjection
+    {
+        public const int FirstZone = 28;
+        public const int LastZone = 38;
+        public const int FirstZoneEpsgCode = 25800 + FirstZone;
+        public const int LastZoneEpsgCode = 25800 + LastZone;
+
+        public int Zone { get; }
+
+        public Etrs89UtmProjection(int zone)
+        {
+            if (zone < FirstZone || zone > LastZone)
+            {
+                throw new ArgumentException($"Invalid ETRS89 UTM zone {zone}.", nameof(zone));
+            }
+
+            Zone = zone;
+            CrsId = $"EPSG:{25800 + Zone}";
+
+            // GRS 1980
+            EquatorialRadius = 6378137d;
+            Flattening = 1d / 298.257222101;
+            ScaleFactor = 0.9996;
+            CentralMeridian = Zone * 6d - 183d;
+            FalseEasting = 5e5;
+            FalseNorthing = 0d;
+        }
+    }
+}

MapControl/Shared/MapProjectionFactory.cs

@@ -35,6 +35,9 @@ namespace MapControl
                 case UpsSouthProjection.DefaultCrsId:
                     return new UpsSouthProjection();
 
+                case Wgs84AutoUtmProjection.DefaultCrsId:
+                    return new Wgs84AutoUtmProjection();
+
                 case OrthographicProjection.DefaultCrsId:
                     return new OrthographicProjection();
 
@@ -50,6 +53,29 @@ namespace MapControl
                 case AzimuthalEquidistantProjection.DefaultCrsId:
                     return new AzimuthalEquidistantProjection();
 
+                default:
+                    return crsId.StartsWith("EPSG:") && int.TryParse(crsId.Substring(5), out int epsgCode)
+                        ? GetProjection(epsgCode)
+                        : null;
+            }
+        }
+
+        public virtual MapProjection GetProjection(int epsgCode)
+        {
+            switch (epsgCode)
+            {
+                case var c when c >= Etrs89UtmProjection.FirstZoneEpsgCode && c <= Etrs89UtmProjection.LastZoneEpsgCode:
+                    return new Etrs89UtmProjection(epsgCode % 100);
+
+                case var c when c >= Nad83UtmProjection.FirstZoneEpsgCode && c <= Nad83UtmProjection.LastZoneEpsgCode:
+                    return new Nad83UtmProjection(epsgCode % 100);
+
+                case var c when c >= Wgs84UtmProjection.FirstZoneNorthEpsgCode && c <= Wgs84UtmProjection.LastZoneNorthEpsgCode:
+                    return new Wgs84UtmProjection(epsgCode % 100, true);
+
+                case var c when c >= Wgs84UtmProjection.FirstZoneSouthEpsgCode && c <= Wgs84UtmProjection.LastZoneSouthEpsgCode:
+                    return new Wgs84UtmProjection(epsgCode % 100, false);
+
                 default:
                     return null;
             }

MapControl/Shared/Nad83UtmProjection.cs

@@ -0,0 +1,40 @@
+// XAML Map Control - https://github.com/ClemensFischer/XAML-Map-Control
+// Copyright © 2024 Clemens Fischer
+// Licensed under the Microsoft Public License (Ms-PL)
+
+using System;
+
+namespace MapControl
+{
+    /// <summary>
+    /// NAD83 UTM Projection with zone number.
+    /// </summary>
+    public class Nad83UtmProjection : TransverseMercatorProjection
+    {
+        public const int FirstZone = 1;
+        public const int LastZone = 23;
+        public const int FirstZoneEpsgCode = 26900 + FirstZone;
+        public const int LastZoneEpsgCode = 26900 + LastZone;
+
+        public int Zone { get; }
+
+        public Nad83UtmProjection(int zone)
+        {
+            if (zone < FirstZone || zone > LastZone)
+            {
+                throw new ArgumentException($"Invalid NAD83 UTM zone {zone}.", nameof(zone));
+            }
+
+            Zone = zone;
+            CrsId = $"EPSG:{26900 + Zone}";
+
+            // GRS 1980
+            EquatorialRadius = 6378137d;
+            Flattening = 1d / 298.257222101;
+            ScaleFactor = 0.9996;
+            CentralMeridian = Zone * 6d - 183d;
+            FalseEasting = 5e5;
+            FalseNorthing = 0d;
+        }
+    }
+}

MapControl/Shared/TransverseMercatorProjection.cs

@@ -0,0 +1,130 @@
+// XAML Map Control - https://github.com/ClemensFischer/XAML-Map-Control
+// Copyright © 2024 Clemens Fischer
+// Licensed under the Microsoft Public License (Ms-PL)
+
+using System;
+#if WPF
+using System.Windows;
+#endif
+
+namespace MapControl
+{
+    /// <summary>
+    /// Universal Transverse Mercator Projection.
+    /// See https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system.
+    /// </summary>
+    public class TransverseMercatorProjection : MapProjection
+    {
+#if NETFRAMEWORK || UWP
+        private static double Math_Atanh(double x) => Math.Log((1d + x) / (1d - x)) / 2d;
+#else
+        private static double Math_Atanh(double x) => Math.Atanh(x);
+#endif
+        public double EquatorialRadius { get; set; } = Wgs84EquatorialRadius;
+        public double Flattening { get; set; } = Wgs84Flattening;
+        public double ScaleFactor { get; set; } = 0.9996;
+        public double CentralMeridian { get; set; }
+        public double FalseEasting { get; set; } = 5e5;
+        public double FalseNorthing { get; set; }
+
+        public TransverseMercatorProjection()
+        {
+            Type = MapProjectionType.TransverseCylindrical;
+        }
+
+        public override Point GetRelativeScale(Location location)
+        {
+            // Precise enough, avoid calculations as in LocationToMap.
+            //
+            return new Point(ScaleFactor, ScaleFactor);
+        }
+
+        public override Point? LocationToMap(Location location)
+        {
+            var n = Flattening / (2d - Flattening);
+            var n2 = n * n;
+            var n3 = n * n2;
+
+            var alpha1 = n / 2d - n2 * 2d / 3d + n3 * 5d / 16d;
+            var alpha2 = n2 * 13d / 48d - n3 * 3d / 5d;
+            var alpha3 = n3 * 61d / 240d;
+
+            var k0A = ScaleFactor * EquatorialRadius / (1d + n) * (1d + n2 / 4d + n2 * n2 / 64d);
+
+            // φ
+            var phi = location.Latitude * Math.PI / 180d;
+
+            // (λ - λ0)
+            var lambda = (location.Longitude - CentralMeridian) * Math.PI / 180d;
+
+            var s = 2d * Math.Sqrt(n) / (1d + n);
+            var sinLat = Math.Sin(phi);
+
+            var t = Math.Sinh(Math_Atanh(sinLat) - s * Math_Atanh(s * sinLat));
+
+            // ξ'
+            var xi_ = Math.Atan(t / Math.Cos(lambda));
+
+            // η'
+            var eta_ = Math_Atanh(Math.Sin(lambda) / Math.Sqrt(1d + t * t));
+
+            var N = FalseNorthing + k0A * (xi_
+                + alpha1 * Math.Sin(2d * xi_) * Math.Cosh(2d * eta_)
+                + alpha2 * Math.Sin(4d * xi_) * Math.Cosh(4d * eta_)
+                + alpha3 * Math.Sin(6d * xi_) * Math.Cosh(6d * eta_));
+
+            var E = FalseEasting + k0A * (eta_
+                + alpha1 * Math.Cos(2d * xi_) * Math.Sinh(2d * eta_)
+                + alpha2 * Math.Cos(4d * xi_) * Math.Sinh(4d * eta_)
+                + alpha3 * Math.Cos(6d * xi_) * Math.Sinh(6d * eta_));
+
+            return new Point(E, N);
+        }
+
+        public override Location MapToLocation(Point point)
+        {
+            var n = Flattening / (2d - Flattening);
+            var n2 = n * n;
+            var n3 = n * n2;
+
+            var beta1 = n / 2d - n2 * 2d / 3d + n3 * 37d / 96d;
+            var beta2 = n2 / 48d + n3 / 15d;
+            var beta3 = n3 * 17d / 480d;
+
+            var delta1 = n * 2d - n2 * 2d / 3d - n3 * 2d;
+            var delta2 = n2 * 7d / 3d - n3 * 8d / 5d;
+            var delta3 = n3 * 56d / 15d;
+
+            var k0A = ScaleFactor * EquatorialRadius / (1d + n) * (1d + n2 / 4d + n2 * n2 / 64d);
+
+            // ξ
+            var xi = (point.Y - FalseNorthing) / k0A;
+
+            // η
+            var eta = (point.X - FalseEasting) / k0A;
+
+            // ξ'
+            var xi_ = xi
+                - beta1 * Math.Sin(2d * xi) * Math.Cosh(2d * eta)
+                - beta2 * Math.Sin(4d * xi) * Math.Cosh(4d * eta)
+                - beta3 * Math.Sin(6d * xi) * Math.Cosh(6d * eta);
+
+            // η'
+            var eta_ = eta
+                - beta1 * Math.Cos(2d * xi) * Math.Sinh(2d * eta)
+                - beta2 * Math.Cos(4d * xi) * Math.Sinh(4d * eta)
+                - beta3 * Math.Cos(6d * xi) * Math.Sinh(6d * eta);
+
+            var chi = Math.Asin(Math.Sin(xi_) / Math.Cosh(eta_));
+
+            var phi = chi
+                + delta1 * Math.Sin(2d * chi)
+                + delta2 * Math.Sin(4d * chi)
+                + delta3 * Math.Sin(6d * chi);
+
+            var lambda = Math.Atan(Math.Sinh(eta_) / Math.Cos(xi_));
+
+            return new Location(phi * 180d / Math.PI, lambda * 180d / Math.PI + CentralMeridian);
+        }
+    }
+}

MapControl/Shared/Wgs84UtmProjection.cs

@@ -0,0 +1,94 @@
+// XAML Map Control - https://github.com/ClemensFischer/XAML-Map-Control
+// Copyright © 2024 Clemens Fischer
+// Licensed under the Microsoft Public License (Ms-PL)
+
+using System;
+
+namespace MapControl
+{
+    /// <summary>
+    /// WGS84 UTM Projection with zone number and north/south flag.
+    /// </summary>
+    public class Wgs84UtmProjection : TransverseMercatorProjection
+    {
+        public const int FirstZone = 1;
+        public const int LastZone = 60;
+        public const int FirstZoneNorthEpsgCode = 32600 + FirstZone;
+        public const int LastZoneNorthEpsgCode = 32600 + LastZone;
+        public const int FirstZoneSouthEpsgCode = 32700 + FirstZone;
+        public const int LastZoneSouthEpsgCode = 32700 + LastZone;
+
+        public int Zone { get; private set; }
+        public bool IsNorth { get; private set; }
+
+        public Wgs84UtmProjection(int zone, bool north)
+        {
+            SetZone(zone, north);
+        }
+
+        protected void SetZone(int zone, bool north)
+        {
+            if (zone < FirstZone || zone > LastZone)
+            {
+                throw new ArgumentException($"Invalid WGS84 UTM zone {zone}.", nameof(zone));
+            }
+
+            Zone = zone;
+            IsNorth = north;
+            CrsId = $"EPSG:{(north ? 32600 : 32700) + zone}";
+            EquatorialRadius = Wgs84EquatorialRadius;
+            Flattening = Wgs84Flattening;
+            ScaleFactor = 0.9996;
+            CentralMeridian = Zone * 6d - 183d;
+            FalseEasting = 5e5;
+            FalseNorthing = IsNorth ? 0d : 1e7;
+        }
+    }
+
+    /// <summary>
+    /// WGS84 UTM Projection with automatic zone selection from projection center.
+    /// </summary>
+    public class Wgs84AutoUtmProjection : Wgs84UtmProjection
+    {
+        public const string DefaultCrsId = "AUTO2:42001";
+
+        private readonly string autoCrsId;
+
+        public Wgs84AutoUtmProjection(string crsId = DefaultCrsId)
+            : base(31, true)
+        {
+            autoCrsId = crsId;
+
+            if (!string.IsNullOrEmpty(autoCrsId))
+            {
+                CrsId = autoCrsId;
+            }
+        }
+
+        public override Location Center
+        {
+            get => base.Center;
+            protected internal set
+            {
+                if (!Equals(base.Center, value))
+                {
+                    base.Center = value;
+
+                    var lon = Location.NormalizeLongitude(value.Longitude);
+                    var zone = (int)Math.Floor(lon / 6d) + 31;
+                    var north = value.Latitude >= 0d;
+
+                    if (Zone != zone || IsNorth != north)
+                    {
+                        SetZone(zone, north);
+
+                        if (!string.IsNullOrEmpty(autoCrsId))
+                        {
+                            CrsId = autoCrsId;
+                        }
+                    }
+                }
+            }
+        }
+    }
+}

MapControl/UWP/MapControl.UWP.csproj

@@ -68,6 +68,9 @@
     <Compile Include="..\Shared\EquirectangularProjection.cs">
       <Link>EquirectangularProjection.cs</Link>
     </Compile>
+    <Compile Include="..\Shared\Etrs89UtmProjection.cs">
+      <Link>Etrs89UtmProjection.cs</Link>
+    </Compile>
     <Compile Include="..\Shared\FilePath.cs">
       <Link>FilePath.cs</Link>
     </Compile>
@@ -140,6 +143,9 @@
     <Compile Include="..\Shared\MapTileLayerBase.cs">
       <Link>MapTileLayerBase.cs</Link>
     </Compile>
+    <Compile Include="..\Shared\Nad83UtmProjection.cs">
+      <Link>Nad83UtmProjection.cs</Link>
+    </Compile>
     <Compile Include="..\Shared\OrthographicProjection.cs">
       <Link>OrthographicProjection.cs</Link>
     </Compile>
@@ -173,6 +179,9 @@
     <Compile Include="..\Shared\Timer.cs">
       <Link>Timer.cs</Link>
     </Compile>
+    <Compile Include="..\Shared\TransverseMercatorProjection.cs">
+      <Link>TransverseMercatorProjection.cs</Link>
+    </Compile>
     <Compile Include="..\Shared\ViewportChangedEventArgs.cs">
       <Link>ViewportChangedEventArgs.cs</Link>
     </Compile>
@@ -182,6 +191,9 @@
     <Compile Include="..\Shared\WebMercatorProjection.cs">
       <Link>WebMercatorProjection.cs</Link>
     </Compile>
+    <Compile Include="..\Shared\Wgs84UtmProjection.cs">
+      <Link>Wgs84UtmProjection.cs</Link>
+    </Compile>
     <Compile Include="..\Shared\WmsImageLayer.cs">
       <Link>WmsImageLayer.cs</Link>
     </Compile>

MapProjections/Shared/GeoApiProjection.cs

@@ -100,6 +100,13 @@ namespace MapControl.Projections
 
         public IMathTransform MapToLocationTransform { get; private set; }
 
+        public override Point GetRelativeScale(Location location)
+        {
+            var k = coordinateSystem?.Projection?.GetParameter("scale_factor")?.Value ?? 1d;
+
+            return new Point(k, k);
+        }
+
         public override Point? LocationToMap(Location location)
         {
             if (LocationToMapTransform == null)

MapProjections/Shared/GeoApiProjectionFactory.cs

@@ -21,36 +21,23 @@ namespace MapControl.Projections
 
         public override MapProjection GetProjection(string crsId)
         {
-            MapProjection projection = null;
-
             switch (crsId)
             {
                 case MapControl.WebMercatorProjection.DefaultCrsId:
-                    projection = new WebMercatorProjection();
-                    break;
+                    return new WebMercatorProjection();
 
                 case MapControl.WorldMercatorProjection.DefaultCrsId:
-                    projection = new WorldMercatorProjection();
-                    break;
+                    return new WorldMercatorProjection();
 
-                case Wgs84AutoUtmProjection.DefaultCrsId:
-                    projection = new Wgs84AutoUtmProjection();
-                    break;
+                case MapControl.Wgs84AutoUtmProjection.DefaultCrsId:
+                    return new Wgs84AutoUtmProjection();
 
                 default:
-                    projection = base.GetProjection(crsId);
-                    break;
+                    return base.GetProjection(crsId);
             }
-
-            if (projection == null && crsId.StartsWith("EPSG:") && int.TryParse(crsId.Substring(5), out int epsgCode))
-            {
-                projection = GetProjection(epsgCode);
-            }
-
-            return projection;
         }
 
-        public virtual MapProjection GetProjection(int epsgCode)
+        public override MapProjection GetProjection(int epsgCode)
         {
             switch (epsgCode)
             {
@@ -75,7 +62,7 @@ namespace MapControl.Projections
                 default:
                     return CoordinateSystemWkts.TryGetValue(epsgCode, out string wkt)
                         ? new GeoApiProjection(wkt)
-                        : null;
+                        : base.GetProjection(epsgCode);
             }
         }
 

MapProjections/Shared/Wgs84UtmProjection.cs

@@ -45,11 +45,9 @@ namespace MapControl.Projections
     /// </summary>
     public class Wgs84AutoUtmProjection : Wgs84UtmProjection
     {
-        public const string DefaultCrsId = "AUTO2:42001";
-
         private readonly string autoCrsId;
 
-        public Wgs84AutoUtmProjection(string crsId = DefaultCrsId)
+        public Wgs84AutoUtmProjection(string crsId = MapControl.Wgs84AutoUtmProjection.DefaultCrsId)
             : base(31, true)
         {
             autoCrsId = crsId;