diff --git a/MapControl/Shared/Etrs89UtmProjection.cs b/MapControl/Shared/Etrs89UtmProjection.cs
new file mode 100644
index 00000000..d8c1ce3b
--- /dev/null
+++ b/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
+{
+ ///
+ /// ETRS89 UTM Projection with zone number.
+ ///
+ 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;
+ }
+ }
+}
diff --git a/MapControl/Shared/MapProjectionFactory.cs b/MapControl/Shared/MapProjectionFactory.cs
index 654de600..b0fe885d 100644
--- a/MapControl/Shared/MapProjectionFactory.cs
+++ b/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;
}
diff --git a/MapControl/Shared/Nad83UtmProjection.cs b/MapControl/Shared/Nad83UtmProjection.cs
new file mode 100644
index 00000000..baad67bd
--- /dev/null
+++ b/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
+{
+ ///
+ /// NAD83 UTM Projection with zone number.
+ ///
+ 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;
+ }
+ }
+}
diff --git a/MapControl/Shared/TransverseMercatorProjection.cs b/MapControl/Shared/TransverseMercatorProjection.cs
new file mode 100644
index 00000000..663a41ce
--- /dev/null
+++ b/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
+{
+ ///
+ /// Universal Transverse Mercator Projection.
+ /// See https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system.
+ ///
+ 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);
+ }
+ }
+}
diff --git a/MapControl/Shared/Wgs84UtmProjection.cs b/MapControl/Shared/Wgs84UtmProjection.cs
new file mode 100644
index 00000000..76014350
--- /dev/null
+++ b/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
+{
+ ///
+ /// WGS84 UTM Projection with zone number and north/south flag.
+ ///
+ 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;
+ }
+ }
+
+ ///
+ /// WGS84 UTM Projection with automatic zone selection from projection center.
+ ///
+ 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;
+ }
+ }
+ }
+ }
+ }
+ }
+}
diff --git a/MapControl/UWP/MapControl.UWP.csproj b/MapControl/UWP/MapControl.UWP.csproj
index 703f653d..2da96286 100644
--- a/MapControl/UWP/MapControl.UWP.csproj
+++ b/MapControl/UWP/MapControl.UWP.csproj
@@ -68,6 +68,9 @@
EquirectangularProjection.cs
+
+ Etrs89UtmProjection.cs
+
FilePath.cs
@@ -140,6 +143,9 @@
MapTileLayerBase.cs
+
+ Nad83UtmProjection.cs
+
OrthographicProjection.cs
@@ -173,6 +179,9 @@
Timer.cs
+
+ TransverseMercatorProjection.cs
+
ViewportChangedEventArgs.cs
@@ -182,6 +191,9 @@
WebMercatorProjection.cs
+
+ Wgs84UtmProjection.cs
+
WmsImageLayer.cs
diff --git a/MapProjections/Shared/GeoApiProjection.cs b/MapProjections/Shared/GeoApiProjection.cs
index ba66e1ff..856d5bbd 100644
--- a/MapProjections/Shared/GeoApiProjection.cs
+++ b/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)
diff --git a/MapProjections/Shared/GeoApiProjectionFactory.cs b/MapProjections/Shared/GeoApiProjectionFactory.cs
index 2e465b61..2c327ad2 100644
--- a/MapProjections/Shared/GeoApiProjectionFactory.cs
+++ b/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);
}
}
diff --git a/MapProjections/Shared/Wgs84UtmProjection.cs b/MapProjections/Shared/Wgs84UtmProjection.cs
index df7fd64b..8360e5b4 100644
--- a/MapProjections/Shared/Wgs84UtmProjection.cs
+++ b/MapProjections/Shared/Wgs84UtmProjection.cs
@@ -45,11 +45,9 @@ namespace MapControl.Projections
///
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;