From e45180b26a869417021f8350a03192a68f241602 Mon Sep 17 00:00:00 2001 From: ClemensFischer Date: Tue, 13 Dec 2022 18:22:18 +0100 Subject: [PATCH] Add map projections Add TransverseMercatorProjection, PolarStereographicProjection and derived UTM/UPS projections to the base library. --- .../Shared/EquirectangularProjection.cs | 3 +- MapControl/Shared/MapProjection.cs | 13 +- MapControl/Shared/MapProjectionFactory.cs | 64 +++++-- .../Shared/PolarStereographicProjection.cs | 144 +++++++++++++++ .../Shared/TransverseMercatorProjection.cs | 169 ++++++++++++++++++ MapControl/Shared/UtmProjection.cs | 132 ++++++++++++++ MapControl/Shared/WebMercatorProjection.cs | 3 +- MapControl/Shared/WorldMercatorProjection.cs | 36 ++-- MapControl/UWP/MapControl.UWP.csproj | 9 + MapProjections/Shared/Ed50UtmProjection.cs | 10 +- MapProjections/Shared/Etrs89UtmProjection.cs | 10 +- MapProjections/Shared/GeoApiProjection.cs | 60 ++----- .../Shared/GeoApiProjectionFactory.cs | 124 +++++++------ .../Shared/PolarStereographicProjection.cs | 142 --------------- .../Shared/WebMercatorProjection.cs | 3 - ...rojection.cs => Wgs84AutoUtmProjection.cs} | 45 +++-- .../Shared/WorldMercatorProjection.cs | 3 - MapProjections/UWP/MapProjections.UWP.csproj | 9 +- 18 files changed, 647 insertions(+), 332 deletions(-) create mode 100644 MapControl/Shared/PolarStereographicProjection.cs create mode 100644 MapControl/Shared/TransverseMercatorProjection.cs create mode 100644 MapControl/Shared/UtmProjection.cs delete mode 100644 MapProjections/Shared/PolarStereographicProjection.cs rename MapProjections/Shared/{AutoUtmProjection.cs => Wgs84AutoUtmProjection.cs} (53%) diff --git a/MapControl/Shared/EquirectangularProjection.cs b/MapControl/Shared/EquirectangularProjection.cs index 2c4127bd..c22edb8d 100644 --- a/MapControl/Shared/EquirectangularProjection.cs +++ b/MapControl/Shared/EquirectangularProjection.cs @@ -17,7 +17,8 @@ namespace MapControl /// public class EquirectangularProjection : MapProjection { - public const string DefaultCrsId = "EPSG:4326"; + public const int DefaultEpsgCode = 4326; + public static readonly string DefaultCrsId = $"EPSG:{DefaultEpsgCode}"; public EquirectangularProjection() { diff --git a/MapControl/Shared/MapProjection.cs b/MapControl/Shared/MapProjection.cs index 41014c28..fb5330c8 100644 --- a/MapControl/Shared/MapProjection.cs +++ b/MapControl/Shared/MapProjection.cs @@ -44,15 +44,12 @@ namespace MapControl /// /// Gets or sets an optional projection center. /// - public Location Center { get; set; } = new Location(); + public virtual Location Center { get; set; } = new Location(); /// /// Gets the relative map scale at the specified Location. /// - public virtual Scale GetRelativeScale(Location location) - { - return new Scale(1d, 1d); - } + public virtual Scale GetRelativeScale(Location location) => new Scale(1d, 1d); /// /// Transforms a Location in geographic coordinates to a Point in projected map coordinates. @@ -137,11 +134,15 @@ namespace MapControl // const double p = 0.01; - return string.Format(CultureInfo.InvariantCulture, "{0},{1},{2},{3}", + var bbox = string.Format(CultureInfo.InvariantCulture, "{0:F2},{1:F2},{2:F2},{3:F2}", p * Math.Ceiling(mapRect.XMin / p), p * Math.Ceiling(mapRect.YMin / p), p * Math.Floor(mapRect.XMax / p), p * Math.Floor(mapRect.YMax / p)); + + System.Diagnostics.Debug.WriteLine(bbox); + + return bbox; } } } diff --git a/MapControl/Shared/MapProjectionFactory.cs b/MapControl/Shared/MapProjectionFactory.cs index b111a030..e2432e4e 100644 --- a/MapControl/Shared/MapProjectionFactory.cs +++ b/MapControl/Shared/MapProjectionFactory.cs @@ -6,22 +6,64 @@ namespace MapControl { public class MapProjectionFactory { + public virtual MapProjection GetProjection(int epsgCode) + { + MapProjection projection = null; + + switch (epsgCode) + { + case WorldMercatorProjection.DefaultEpsgCode: + projection = new WorldMercatorProjection(); + break; + + case WebMercatorProjection.DefaultEpsgCode: + projection = new WebMercatorProjection(); + break; + + case EquirectangularProjection.DefaultEpsgCode: + projection = new EquirectangularProjection(); + break; + + case UpsNorthProjection.DefaultEpsgCode: + projection = new UpsNorthProjection(); + break; + + case UpsSouthProjection.DefaultEpsgCode: + projection = new UpsSouthProjection(); + break; + + case var c when c >= Etrs89UtmProjection.FirstZoneEpsgCode && c <= Etrs89UtmProjection.LastZoneEpsgCode: + projection = new Etrs89UtmProjection(epsgCode % 100); + break; + + case var c when c >= Wgs84UtmProjection.FirstZoneNorthEpsgCode && c <= Wgs84UtmProjection.LastZoneNorthEpsgCode: + projection = new Wgs84UtmProjection(epsgCode % 100, true); + break; + + case var c when c >= Wgs84UtmProjection.FirstZoneSouthEpsgCode && c <= Wgs84UtmProjection.LastZoneSouthEpsgCode: + projection = new Wgs84UtmProjection(epsgCode % 100, false); + break; + + default: + break; + } + + return projection; + } + public virtual MapProjection GetProjection(string crsId) { + if (crsId.StartsWith("EPSG:") && int.TryParse(crsId.Substring(5), out int epsgCode)) + { + return GetProjection(epsgCode); + } + MapProjection projection = null; switch (crsId) { - case WorldMercatorProjection.DefaultCrsId: - projection = new WorldMercatorProjection(); - break; - - case WebMercatorProjection.DefaultCrsId: - projection = new WebMercatorProjection(); - break; - - case EquirectangularProjection.DefaultCrsId: - projection = new EquirectangularProjection(); + case Wgs84AutoUtmProjection.DefaultCrsId: + projection = new Wgs84AutoUtmProjection(); break; case OrthographicProjection.DefaultCrsId: @@ -40,7 +82,7 @@ namespace MapControl projection = new StereographicProjection(); break; - case "EPSG:97003": // proprietary CRS ID + case "AUTO2:97003": // proprietary CRS ID projection = new AzimuthalEquidistantProjection(crsId); break; diff --git a/MapControl/Shared/PolarStereographicProjection.cs b/MapControl/Shared/PolarStereographicProjection.cs new file mode 100644 index 00000000..4f5b0871 --- /dev/null +++ b/MapControl/Shared/PolarStereographicProjection.cs @@ -0,0 +1,144 @@ +// XAML Map Control - https://github.com/ClemensFischer/XAML-Map-Control +// © 2022 Clemens Fischer +// Licensed under the Microsoft Public License (Ms-PL) + +using System; +#if !WINUI && !UWP +using System.Windows; +#endif + +namespace MapControl +{ + /// + /// Elliptical Polar Stereographic Projection with a given scale factor at the pole and + /// optional false easting and northing, as used by the UPS North and UPS South projections. + /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/pp/1395/report.pdf), p.154-163. + /// + public class PolarStereographicProjection : MapProjection + { + public PolarStereographicProjection() + { + Type = MapProjectionType.Azimuthal; + } + + public double EquatorialRadius { get; set; } = Wgs84EquatorialRadius; + public double Flattening { get; set; } = Wgs84Flattening; + public double ScaleFactor { get; set; } = 0.994; + public double FalseEasting { get; set; } = 2e6; + public double FalseNorthing { get; set; } = 2e6; + public bool IsNorth { get; set; } + + public override Scale GetRelativeScale(Location location) + { + var lat = location.Latitude * Math.PI / 180d; + + if (!IsNorth) + { + lat = -lat; + } + + var e = Math.Sqrt((2d - Flattening) * Flattening); + var eSinLat = e * Math.Sin(lat); + + var t = Math.Tan(Math.PI / 4d - lat / 2d) + / Math.Pow((1d - eSinLat) / (1d + eSinLat), e / 2d); // p.161 (15-9) + var r = 2d * EquatorialRadius * ScaleFactor * t + / Math.Sqrt(Math.Pow(1d + e, 1d + e) * Math.Pow(1d - e, 1d - e)); // p.161 (21-33) + + var m = Math.Cos(lat) / Math.Sqrt(1d - eSinLat * eSinLat); // p.160 (14-15) + var k = r / (EquatorialRadius * m); // p.161 (21-32) + + return new Scale(k, k); + } + + public override Point? LocationToMap(Location location) + { + var lat = location.Latitude * Math.PI / 180d; + var lon = location.Longitude * Math.PI / 180d; + + if (!IsNorth) + { + lat = -lat; + lon = -lon; + } + + var e = Math.Sqrt((2d - Flattening) * Flattening); + var eSinLat = e * Math.Sin(lat); + + var t = Math.Tan(Math.PI / 4d - lat / 2d) + / Math.Pow((1d - eSinLat) / (1d + eSinLat), e / 2d); // p.161 (15-9) + var r = 2d * EquatorialRadius * ScaleFactor * t + / Math.Sqrt(Math.Pow(1d + e, 1d + e) * Math.Pow(1d - e, 1d - e)); // p.161 (21-33) + + var x = r * Math.Sin(lon); // p.161 (21-30) + var y = -r * Math.Cos(lon); // p.161 (21-31) + + if (!IsNorth) + { + x = -x; + y = -y; + } + + return new Point(x + FalseEasting, y + FalseNorthing); + } + + public override Location MapToLocation(Point point) + { + var x = point.X - FalseEasting; + var y = point.Y - FalseNorthing; + + if (!IsNorth) + { + x = -x; + y = -y; + } + + var e = Math.Sqrt((2d - Flattening) * Flattening); + + 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 lat = WorldMercatorProjection.LatitudeFromSeriesApproximation(e, t); // p.162 (3-5) + var lon = Math.Atan2(x, -y); // p.162 (20-16) + + if (!IsNorth) + { + lat = -lat; + lon = -lon; + } + + return new Location(lat * 180d / Math.PI, lon * 180d / Math.PI); + } + } + + /// + /// Universal Polar Stereographic North Projection - EPSG:32661. + /// + public class UpsNorthProjection : PolarStereographicProjection + { + public const int DefaultEpsgCode = 32661; + public static readonly string DefaultCrsId = $"EPSG:{DefaultEpsgCode}"; + + public UpsNorthProjection() + { + CrsId = DefaultCrsId; + IsNorth = true; + } + } + + /// + /// Universal Polar Stereographic South Projection - EPSG:32761. + /// + public class UpsSouthProjection : PolarStereographicProjection + { + public const int DefaultEpsgCode = 32761; + public static readonly string DefaultCrsId = $"EPSG:{DefaultEpsgCode}"; + + public UpsSouthProjection() + { + CrsId = DefaultCrsId; + IsNorth = false; + } + } +} diff --git a/MapControl/Shared/TransverseMercatorProjection.cs b/MapControl/Shared/TransverseMercatorProjection.cs new file mode 100644 index 00000000..da7cf989 --- /dev/null +++ b/MapControl/Shared/TransverseMercatorProjection.cs @@ -0,0 +1,169 @@ +// XAML Map Control - https://github.com/ClemensFischer/XAML-Map-Control +// © 2022 Clemens Fischer +// Licensed under the Microsoft Public License (Ms-PL) + +using System; +#if !WINUI && !UWP +using System.Windows; +#endif + +namespace MapControl +{ + /// + /// Universal Transverse Mercator Projection. + /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/pp/1395/report.pdf), p.57-64. + /// + public class TransverseMercatorProjection : MapProjection + { + public TransverseMercatorProjection() + { + Type = MapProjectionType.TransverseCylindrical; + } + + 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; } + public double FalseNorthing { get; set; } + + public override Scale GetRelativeScale(Location location) + { + var k = ScaleFactor; + + if (location.Latitude > -90d && location.Latitude < 90d) + { + var lat = location.Latitude * Math.PI / 180d; + var lon = (location.Longitude - CentralMeridian) * Math.PI / 180d; + var cosLat = Math.Cos(lat); + var tanLat = Math.Tan(lat); + + var e_2 = (2d - Flattening) * Flattening; + var E_2 = e_2 / (1d - e_2); // p.61 (8-12) + var T = tanLat * tanLat; // p.61 (8-13) + var C = E_2 * cosLat * cosLat; // p.61 (8-14) + var A = lon * cosLat; // p.61 (8-15) + + var T_2 = T * T; + var C_2 = C * C; + var A_2 = A * A; + var A_4 = A_2 * A_2; + var A_6 = A_2 * A_4; + + k = ScaleFactor * (1d + + (1d + C) * A_2 / 2d + + (5d - 4d * T + 42d * C + 13d * C_2 + 28d * E_2) * A_4 / 24d + + (61d - 148d * T + 16d * T_2) * A_6 / 720d); // p.61 (8-11) + } + + return new Scale(k, k); + } + + public override Point? LocationToMap(Location location) + { + double x, y; + + var lat = location.Latitude * Math.PI / 180d; + var e_2 = (2d - Flattening) * Flattening; + var e_4 = e_2 * e_2; + var e_6 = e_2 * e_4; + var M = EquatorialRadius * (1d - e_2 / 4d - 3d * e_4 / 64d - 5d * e_6 / 256d) * lat; // p.61 (3-21) + + if (lat > -Math.PI / 2d && lat < Math.PI / 2d) + { + var lon = (location.Longitude - CentralMeridian) * Math.PI / 180d; + var sinLat = Math.Sin(lat); + var cosLat = Math.Cos(lat); + var tanLat = Math.Tan(lat); + + var E_2 = e_2 / (1d - e_2); // p.61 (8-12) + var T = tanLat * tanLat; // p.61 (8-13) + var C = E_2 * cosLat * cosLat; // p.61 (8-14) + var A = lon * cosLat; // p.61 (8-15) + + var N = EquatorialRadius / Math.Sqrt(1d - e_2 * sinLat * sinLat); // p.61 (4-20) + + M += EquatorialRadius * + (-(3d * e_2 / 8d + 3d * e_4 / 32d + 45d * e_6 / 1024d) * Math.Sin(2d * lat) + + (15d * e_4 / 256d + 45d * e_6 / 1024d) * Math.Sin(4d * lat) + - (35d * e_6 / 3072d) * Math.Sin(6d * lat)); // p.61 (3-21) + + var T_2 = T * T; + var C_2 = C * C; + var A_2 = A * A; + var A_3 = A * A_2; + var A_4 = A * A_3; + var A_5 = A * A_4; + var A_6 = A * A_5; + + x = ScaleFactor * N * (A + + (1d - T + C) * A_3 / 6d + + (5d - 18d * T + T_2 + 72d * C - 58d * E_2) * A_5 / 120d); // p.61 (8-9) + + y = ScaleFactor * (M + N * tanLat * (A_2 / 2d + + (5d - T + 9d * C + 4d * C_2) * A_4 / 24d + + (61d - 58d * T + T_2 + 600d * C - 330d * E_2) * A_6 / 720d)); // p.61 (8-10) + } + else + { + x = 0d; + y = ScaleFactor * M; + } + + return new Point(x + FalseEasting, y + FalseNorthing); + } + + public override Location MapToLocation(Point point) + { + var x = point.X - FalseEasting; + var y = point.Y - FalseNorthing; + + var e_2 = (2d - Flattening) * Flattening; + var e_4 = e_2 * e_2; + var e_6 = e_2 * e_4; + + var M = y / ScaleFactor; // p.63 (8-20) + var mu = M / (EquatorialRadius * (1d - e_2 / 4d - 3d * e_4 / 64d - 5d * e_6 / 256d)); // p.63 (7-19) + var e1 = (1d - Math.Sqrt(1d - e_2)) / (1d + Math.Sqrt(1d - e_2)); // p.63 (3-24) + var e1_2 = e1 * e1; + var e1_3 = e1 * e1_2; + var e1_4 = e1 * e1_3; + + var lat1 = mu + + (3d * e1 / 2d - 27d * e1_3 / 32d) * Math.Sin(2d * mu) + + (21d * e1_2 / 16d - 55d * e1_4 / 32d) * Math.Sin(4d * mu) + + (151d * e1_3 / 96d) * Math.Sin(6d * mu) + + (1097d * e1_4 / 512d) * Math.Sin(8d * mu); // p.63 (3-26) + + var sinLat1 = Math.Sin(lat1); + var cosLat1 = Math.Cos(lat1); + var tanLat1 = Math.Tan(lat1); + + var E_2 = e_2 / (1d - e_2); // p.64 (8-12) + var C1 = E_2 * cosLat1 * cosLat1; // p.64 (8-21) + var T1 = tanLat1 * tanLat1; // p.64 (8-22) + var N1 = EquatorialRadius / Math.Sqrt(1d - e_2 * sinLat1 * sinLat1); // p.64 (8-23) + var R1 = EquatorialRadius * (1d - e_2) / Math.Pow(1d - e_2 * sinLat1 * sinLat1, 1.5); // p.64 (8-24) + var D = x / (N1 * ScaleFactor); // p.64 (8-25) + + var C1_2 = C1 * C1; + var T1_2 = T1 * T1; + var D_2 = D * D; + var D_3 = D * D_2; + var D_4 = D * D_3; + var D_5 = D * D_4; + var D_6 = D * D_5; + + var lat = lat1 - (N1 * tanLat1 / R1) * (D_2 / 2d + - (5d + 3d * T1 + 10d * C1 - 4d * C1_2 - 9d * E_2) * D_4 / 24d + + (61d + 90d * T1 + 298d * C1 + 45d * T1_2 - 252d * E_2 - 3d * C1_2) * D_6 / 720d); // p.63 (8-17) + + var lon = (D + - (1d + 2d * T1 + C1) * D_3 / 6d + + (5d - 2d * C1 + 28d * T1 - 3d * C1_2 + 8d * E_2 + 24d * T1_2) * D_5 / 120d) + / cosLat1; // p.63 (8-18) + + return new Location(lat * 180d / Math.PI, lon * 180d / Math.PI + CentralMeridian); + } + } +} diff --git a/MapControl/Shared/UtmProjection.cs b/MapControl/Shared/UtmProjection.cs new file mode 100644 index 00000000..90daba25 --- /dev/null +++ b/MapControl/Shared/UtmProjection.cs @@ -0,0 +1,132 @@ +// XAML Map Control - https://github.com/ClemensFischer/XAML-Map-Control +// © 2022 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; + + private int zone; + + public Etrs89UtmProjection(int zone) + { + Zone = zone; + } + + public int Zone + { + get => zone; + set + { + if (value < FirstZone || value > LastZone) + { + throw new ArgumentException($"Invalid ETRS89 UTM zone {value}.", nameof(value)); + } + + zone = value; + CrsId = $"EPSG:{zone - FirstZone + FirstZoneEpsgCode}"; + EquatorialRadius = 6378137d; + Flattening = 1 / 298.257222101; // GRS 1980 + ScaleFactor = 0.9996; + CentralMeridian = zone * 6d - 183d; + FalseEasting = 5e5; + FalseNorthing = 0d; + } + } + } + + /// + /// 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; + + protected Wgs84UtmProjection() + { + } + + public Wgs84UtmProjection(int zone, bool north) + { + SetZone(zone, north); + } + + public int Zone { get; private set; } + public bool IsNorth { get; private set; } + + public void SetZone(int zone, bool north, string crsId = null) + { + if (zone < FirstZone || zone > LastZone) + { + throw new ArgumentException($"Invalid WGS84 UTM zone {zone}.", nameof(zone)); + } + + var epsgCode = zone - FirstZone + (north ? FirstZoneNorthEpsgCode : FirstZoneSouthEpsgCode); + + Zone = zone; + IsNorth = north; + CrsId = crsId ?? $"EPSG:{epsgCode}"; + 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"; + + public Wgs84AutoUtmProjection(bool useZoneCrsId = false) + { + UseZoneCrsId = useZoneCrsId; + UpdateZone(); + } + + public bool UseZoneCrsId { get; } + + public override Location Center + { + get => base.Center; + set + { + if (!Equals(base.Center, value)) + { + base.Center = value; + UpdateZone(); + } + } + } + + private void UpdateZone() + { + var lon = Location.NormalizeLongitude(Center.Longitude); + var zone = (int)Math.Floor(lon / 6d) + 31; + var north = Center.Latitude >= 0d; + + if (Zone != zone || IsNorth != north || string.IsNullOrEmpty(CrsId)) + { + SetZone(zone, north, UseZoneCrsId ? null : DefaultCrsId); + } + } + + } +} diff --git a/MapControl/Shared/WebMercatorProjection.cs b/MapControl/Shared/WebMercatorProjection.cs index 8e123b20..92fc2060 100644 --- a/MapControl/Shared/WebMercatorProjection.cs +++ b/MapControl/Shared/WebMercatorProjection.cs @@ -15,7 +15,8 @@ namespace MapControl /// public class WebMercatorProjection : MapProjection { - public const string DefaultCrsId = "EPSG:3857"; + public const int DefaultEpsgCode = 3857; + public static readonly string DefaultCrsId = $"EPSG:{DefaultEpsgCode}"; public WebMercatorProjection() { diff --git a/MapControl/Shared/WorldMercatorProjection.cs b/MapControl/Shared/WorldMercatorProjection.cs index e61cc4f6..2385aa34 100644 --- a/MapControl/Shared/WorldMercatorProjection.cs +++ b/MapControl/Shared/WorldMercatorProjection.cs @@ -15,10 +15,8 @@ namespace MapControl /// public class WorldMercatorProjection : MapProjection { - private const double ConvergenceTolerance = 1e-6; - private const int MaxIterations = 10; - - public const string DefaultCrsId = "EPSG:3395"; + public const int DefaultEpsgCode = 3395; + public static readonly string DefaultCrsId = $"EPSG:{DefaultEpsgCode}"; public WorldMercatorProjection() { @@ -62,31 +60,33 @@ namespace MapControl } var lat = latitude * Math.PI / 180d; + var eSinLat = Wgs84Eccentricity * Math.Sin(lat); + var f = Math.Pow((1d - eSinLat) / (1d + eSinLat), Wgs84Eccentricity / 2d); - return Math.Log(Math.Tan(lat / 2d + Math.PI / 4d) * ConformalFactor(lat)) * 180d / Math.PI; // p.44 (7-7) + return Math.Log(Math.Tan(lat / 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) - var lat = Math.PI / 2d - 2d * Math.Atan(t); // p.44 (7-11) - var relChange = 1d; - for (var i = 0; i < MaxIterations && relChange > ConvergenceTolerance; i++) - { - var newLat = Math.PI / 2d - 2d * Math.Atan(t * ConformalFactor(lat)); // p.44 (7-9) - relChange = Math.Abs(1d - newLat / lat); - lat = newLat; - } - - return lat * 180d / Math.PI; + return LatitudeFromSeriesApproximation(Wgs84Eccentricity, t) * 180d / Math.PI; } - private static double ConformalFactor(double lat) + internal static double LatitudeFromSeriesApproximation(double e, double t) { - var eSinLat = Wgs84Eccentricity * Math.Sin(lat); + var e_2 = e * e; + var e_4 = e_2 * e_2; + var e_6 = e_2 * e_4; + var e_8 = e_2 * e_6; - return Math.Pow((1d - eSinLat) / (1d + eSinLat), Wgs84Eccentricity / 2d); + var lat = Math.PI / 2d - 2d * Math.Atan(t); // p.45 (7-13) + + return lat + + (e_2 / 2d + 5d * e_4 / 24d + e_6 / 12d + 13d * e_8 / 360d) * Math.Sin(2d * lat) + + (7d * e_4 / 48d + 29d * e_6 / 240d + 811d * e_8 / 11520d) * Math.Sin(4d * lat) + + (7d * e_6 / 120d + 81d * e_8 / 1120d) * Math.Sin(6d * lat) + + (4279d * e_8 / 161280d) * Math.Sin(8d * lat); // p.45 (3-5) } } } diff --git a/MapControl/UWP/MapControl.UWP.csproj b/MapControl/UWP/MapControl.UWP.csproj index 28183fe1..61629a28 100644 --- a/MapControl/UWP/MapControl.UWP.csproj +++ b/MapControl/UWP/MapControl.UWP.csproj @@ -146,6 +146,9 @@ OrthographicProjection.cs + + PolarStereographicProjection.cs + PushpinBorder.cs @@ -170,6 +173,12 @@ TileSource.cs + + TransverseMercatorProjection.cs + + + UtmProjection.cs + ViewportChangedEventArgs.cs diff --git a/MapProjections/Shared/Ed50UtmProjection.cs b/MapProjections/Shared/Ed50UtmProjection.cs index 4a4fac7f..3550656b 100644 --- a/MapProjections/Shared/Ed50UtmProjection.cs +++ b/MapProjections/Shared/Ed50UtmProjection.cs @@ -15,8 +15,8 @@ namespace MapControl.Projections throw new ArgumentException($"Invalid UTM zone {zone}.", nameof(zone)); } - const string wktFormat - = "PROJCS[\"ED50 / UTM zone {0}N\"," + CoordinateSystemWkt + = $"PROJCS[\"ED50 / UTM zone {zone}N\"," + "GEOGCS[\"ED50\"," + "DATUM[\"European_Datum_1950\"," + "SPHEROID[\"International 1924\",6378388,297," @@ -30,7 +30,7 @@ namespace MapControl.Projections + "AUTHORITY[\"EPSG\",\"4230\"]]," + "PROJECTION[\"Transverse_Mercator\"]," + "PARAMETER[\"latitude_of_origin\",0]," - + "PARAMETER[\"central_meridian\",{1}]," + + $"PARAMETER[\"central_meridian\",{6 * zone - 183}]," + "PARAMETER[\"scale_factor\",0.9996]," + "PARAMETER[\"false_easting\",500000]," + "PARAMETER[\"false_northing\",0]," @@ -38,9 +38,7 @@ namespace MapControl.Projections + "AUTHORITY[\"EPSG\",\"9001\"]]," + "AXIS[\"Easting\",EAST]," + "AXIS[\"Northing\",NORTH]," - + "AUTHORITY[\"EPSG\",\"230{0}\"]]"; - - CoordinateSystemWkt = string.Format(wktFormat, zone, 6 * zone - 183); + + $"AUTHORITY[\"EPSG\",\"230{zone}\"]]"; } } } diff --git a/MapProjections/Shared/Etrs89UtmProjection.cs b/MapProjections/Shared/Etrs89UtmProjection.cs index 71aa201c..ce707efb 100644 --- a/MapProjections/Shared/Etrs89UtmProjection.cs +++ b/MapProjections/Shared/Etrs89UtmProjection.cs @@ -15,8 +15,8 @@ namespace MapControl.Projections throw new ArgumentException($"Invalid UTM zone {zone}.", nameof(zone)); } - const string wktFormat - = "PROJCS[\"ETRS89 / UTM zone {0}N\"," + CoordinateSystemWkt + = $"PROJCS[\"ETRS89 / UTM zone {zone}N\"," + "GEOGCS[\"ETRS89\"," + "DATUM[\"European_Terrestrial_Reference_System_1989\"," + "SPHEROID[\"GRS 1980\",6378137,298.257222101," @@ -30,7 +30,7 @@ namespace MapControl.Projections + "AUTHORITY[\"EPSG\",\"4258\"]]," + "PROJECTION[\"Transverse_Mercator\"]," + "PARAMETER[\"latitude_of_origin\",0]," - + "PARAMETER[\"central_meridian\",{1}]," + + $"PARAMETER[\"central_meridian\",{6 * zone - 183}]," + "PARAMETER[\"scale_factor\",0.9996]," + "PARAMETER[\"false_easting\",500000]," + "PARAMETER[\"false_northing\",0]," @@ -38,9 +38,7 @@ namespace MapControl.Projections + "AUTHORITY[\"EPSG\",\"9001\"]]," + "AXIS[\"Easting\",EAST]," + "AXIS[\"Northing\",NORTH]," - + "AUTHORITY[\"EPSG\",\"258{0}\"]]"; - - CoordinateSystemWkt = string.Format(wktFormat, zone, 6 * zone - 183); + + $"AUTHORITY[\"EPSG\",\"258{zone}\"]]"; } } } diff --git a/MapProjections/Shared/GeoApiProjection.cs b/MapProjections/Shared/GeoApiProjection.cs index ac236d35..c965af9c 100644 --- a/MapProjections/Shared/GeoApiProjection.cs +++ b/MapProjections/Shared/GeoApiProjection.cs @@ -9,11 +9,7 @@ using ProjNet.CoordinateSystems; using ProjNet.CoordinateSystems.Transformations; using System; using System.Globalization; -#if WINUI -using Windows.Foundation; -#elif UWP -using Windows.Foundation; -#else +#if !WINUI && !UWP using System.Windows; #endif @@ -24,9 +20,7 @@ namespace MapControl.Projections /// public class GeoApiProjection : MapProjection { - private ICoordinateSystem coordinateSystem; - private double scaleFactor; - private string bboxFormat; + private IProjectedCoordinateSystem coordinateSystem; protected GeoApiProjection() { @@ -45,13 +39,13 @@ namespace MapControl.Projections public string CoordinateSystemWkt { get => CoordinateSystem?.WKT; - protected set => CoordinateSystem = new CoordinateSystemFactory().CreateFromWkt(value); + protected set => CoordinateSystem = new CoordinateSystemFactory().CreateFromWkt(value) as IProjectedCoordinateSystem; } /// /// Gets or sets the ICoordinateSystem of the MapProjection. /// - public ICoordinateSystem CoordinateSystem + public IProjectedCoordinateSystem CoordinateSystem { get => coordinateSystem; protected set @@ -72,42 +66,31 @@ namespace MapControl.Projections ? string.Format("{0}:{1}", coordinateSystem.Authority, coordinateSystem.AuthorityCode) : ""; - var projection = (coordinateSystem as IProjectedCoordinateSystem)?.Projection; - - if (projection != null) + if (CrsId == "EPSG:3857") { + Type = MapProjectionType.WebMercator; + } + else + { + var projection = coordinateSystem.Projection + ?? throw new ArgumentException("CoordinateSystem.Projection must not be null.", nameof(value)); + var centralMeridian = projection.GetParameter("central_meridian") ?? projection.GetParameter("longitude_of_origin"); var centralParallel = projection.GetParameter("central_parallel") ?? projection.GetParameter("latitude_of_origin"); var falseEasting = projection.GetParameter("false_easting"); var falseNorthing = projection.GetParameter("false_northing"); - if (CrsId == "EPSG:3857") - { - Type = MapProjectionType.WebMercator; - } - else if ( - (centralMeridian == null || centralMeridian.Value == 0d) && + if ((centralMeridian == null || centralMeridian.Value == 0d) && (centralParallel == null || centralParallel.Value == 0d) && (falseEasting == null || falseEasting.Value == 0d) && (falseNorthing == null || falseNorthing.Value == 0d)) { Type = MapProjectionType.NormalCylindrical; } - else if ( - projection.Name.StartsWith("UTM") || - projection.Name.StartsWith("Transverse")) + else if (projection.Name.StartsWith("UTM") || projection.Name.StartsWith("Transverse")) { Type = MapProjectionType.TransverseCylindrical; } - - scaleFactor = 1d; - bboxFormat = "{0},{1},{2},{3}"; - } - else - { - Type = MapProjectionType.NormalCylindrical; - scaleFactor = Wgs84MeterPerDegree; - bboxFormat = "{1},{0},{3},{2}"; } } } @@ -123,15 +106,14 @@ namespace MapControl.Projections throw new InvalidOperationException("The CoordinateSystem property is not set."); } - var coordinate = LocationToMapTransform.Transform( - new Coordinate(location.Longitude, location.Latitude)); + var coordinate = LocationToMapTransform.Transform(new Coordinate(location.Longitude, location.Latitude)); if (coordinate == null) { return null; } - return new Point(coordinate.X * scaleFactor, coordinate.Y * scaleFactor); + return new Point(coordinate.X, coordinate.Y); } public override Location MapToLocation(Point point) @@ -141,19 +123,15 @@ namespace MapControl.Projections throw new InvalidOperationException("The CoordinateSystem property is not set."); } - var coordinate = MapToLocationTransform.Transform( - new Coordinate(point.X / scaleFactor, point.Y / scaleFactor)); + var coordinate = MapToLocationTransform.Transform(new Coordinate(point.X, point.Y)); return new Location(coordinate.Y, coordinate.X); } public override string GetBboxValue(MapRect rect) { - return string.Format(CultureInfo.InvariantCulture, bboxFormat, - rect.XMin / scaleFactor, - rect.YMin / scaleFactor, - rect.XMax / scaleFactor, - rect.YMax / scaleFactor); + return string.Format(CultureInfo.InvariantCulture, + "{0},{1},{2},{3}", rect.XMin, rect.YMin, rect.XMax, rect.YMax); } } } diff --git a/MapProjections/Shared/GeoApiProjectionFactory.cs b/MapProjections/Shared/GeoApiProjectionFactory.cs index 839a03ba..b6e56647 100644 --- a/MapProjections/Shared/GeoApiProjectionFactory.cs +++ b/MapProjections/Shared/GeoApiProjectionFactory.cs @@ -8,79 +8,75 @@ namespace MapControl.Projections { public class GeoApiProjectionFactory : MapProjectionFactory { - public const int WorldMercator = 3395; - public const int WebMercator = 3857; - public const int AutoUtm = 42001; - public const int Ed50UtmFirst = 23028; - public const int Ed50UtmLast = 23038; - public const int Etrs89UtmFirst = 25828; - public const int Etrs89UtmLast = 25838; - public const int Wgs84UtmNorthFirst = 32601; - public const int Wgs84UtmNorthLast = 32660; - public const int Wgs84UpsNorth = 32661; - public const int Wgs84UtmSouthFirst = 32701; - public const int Wgs84UtmSouthLast = 32760; - public const int Wgs84UpsSouth = 32761; + private const int WorldMercator = 3395; + private const int WebMercator = 3857; + private const int Ed50UtmFirst = 23028; + private const int Ed50UtmLast = 23038; + private const int Etrs89UtmFirst = 25828; + private const int Etrs89UtmLast = 25838; + private const int Wgs84UtmNorthFirst = 32601; + private const int Wgs84UtmNorthLast = 32660; + private const int Wgs84UtmSouthFirst = 32701; + private const int Wgs84UtmSouthLast = 32760; public Dictionary CoordinateSystemWkts { get; } = new Dictionary(); + public override MapProjection GetProjection(int epsgCode) + { + MapProjection projection = null; + + if (CoordinateSystemWkts.TryGetValue(epsgCode, out string wkt)) + { + projection = new GeoApiProjection(wkt); + } + else + { + switch (epsgCode) + { + case WorldMercator: + projection = new WorldMercatorProjection(); + break; + + case WebMercator: + projection = new WebMercatorProjection(); + break; + + case int c when c >= Ed50UtmFirst && c <= Ed50UtmLast: + projection = new Ed50UtmProjection(epsgCode % 100); + break; + + case int c when c >= Etrs89UtmFirst && c <= Etrs89UtmLast: + projection = new Etrs89UtmProjection(epsgCode % 100); + break; + + case int c when c >= Wgs84UtmNorthFirst && c <= Wgs84UtmNorthLast: + projection = new Wgs84UtmProjection(epsgCode % 100, true); + break; + + case int c when c >= Wgs84UtmSouthFirst && c <= Wgs84UtmSouthLast: + projection = new Wgs84UtmProjection(epsgCode % 100, false); + break; + + default: + break; + } + } + + return projection ?? base.GetProjection(epsgCode); + } + public override MapProjection GetProjection(string crsId) { MapProjection projection = null; - var str = crsId.StartsWith("EPSG:") ? crsId.Substring(5) - : crsId.StartsWith("AUTO2:") ? crsId.Substring(6) - : null; - if (int.TryParse(str, out int code)) + switch (crsId) { - if (CoordinateSystemWkts.TryGetValue(code, out string wkt)) - { - projection = new GeoApiProjection(wkt); - } - else - { - switch (code) - { - case WorldMercator: - projection = new WorldMercatorProjection(); - break; + case Wgs84AutoUtmProjection.DefaultCrsId: + projection = new Wgs84AutoUtmProjection(); + break; - case WebMercator: - projection = new WebMercatorProjection(); - break; - - case AutoUtm: - projection = new AutoUtmProjection(); - break; - - case int c when c >= Ed50UtmFirst && c <= Ed50UtmLast: - projection = new Ed50UtmProjection(code % 100); - break; - - case int c when c >= Etrs89UtmFirst && c <= Etrs89UtmLast: - projection = new Etrs89UtmProjection(code % 100); - break; - - case int c when c >= Wgs84UtmNorthFirst && c <= Wgs84UtmNorthLast: - projection = new Wgs84UtmProjection(code % 100, true); - break; - - case int c when c >= Wgs84UtmSouthFirst && c <= Wgs84UtmSouthLast: - projection = new Wgs84UtmProjection(code % 100, false); - break; - - case Wgs84UpsNorth: - projection = new UpsNorthProjection(); - break; - - case Wgs84UpsSouth: - projection = new UpsSouthProjection(); - break; - - default: - break; - } - } + default: + break; } return projection ?? base.GetProjection(crsId); diff --git a/MapProjections/Shared/PolarStereographicProjection.cs b/MapProjections/Shared/PolarStereographicProjection.cs deleted file mode 100644 index 2aab6858..00000000 --- a/MapProjections/Shared/PolarStereographicProjection.cs +++ /dev/null @@ -1,142 +0,0 @@ -// XAML Map Control - https://github.com/ClemensFischer/XAML-Map-Control -// © 2022 Clemens Fischer -// Licensed under the Microsoft Public License (Ms-PL) - -using System; -#if !UWP -using System.Windows; -#endif - -namespace MapControl.Projections -{ - /// - /// Elliptical Polar Stereographic Projection with a given scale factor at the pole and - /// optional false easting and northing, as used by the UPS North and UPS South projections. - /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/pp/1395/report.pdf), p.160-162. - /// - public class PolarStereographicProjection : MapProjection - { - private const double ConvergenceTolerance = 1e-6; - private const int MaxIterations = 10; - - public bool IsNorth { get; } - public double ScaleFactor { get; } - public double FalseEasting { get; } - public double FalseNorthing { get; } - - public PolarStereographicProjection(bool isNorth, double scaleFactor, double falseEasting, double falseNorthing) - { - Type = MapProjectionType.Azimuthal; - IsNorth = isNorth; - ScaleFactor = scaleFactor; - FalseEasting = falseEasting; - FalseNorthing = falseNorthing; - } - - public override Scale GetRelativeScale(Location location) - { - var lat = (IsNorth ? location.Latitude : -location.Latitude) * Math.PI / 180d; - var a = Wgs84EquatorialRadius; - var e = Wgs84Eccentricity; - var s = Math.Sqrt(Math.Pow(1 + e, 1 + e) * Math.Pow(1 - e, 1 - e)); - var t = Math.Tan(Math.PI / 4d - lat / 2d) / ConformalFactor(lat); - var rho = 2d * a * ScaleFactor * t / s; - var eSinLat = e * Math.Sin(lat); - var m = Math.Cos(lat) / Math.Sqrt(1d - eSinLat * eSinLat); - var k = rho / (a * m); - - return new Scale(k, k); - } - - public override Point? LocationToMap(Location location) - { - var lat = location.Latitude * Math.PI / 180d; - var lon = location.Longitude * Math.PI / 180d; - - if (IsNorth) - { - lon = Math.PI - lon; - } - else - { - lat = -lat; - } - - var a = Wgs84EquatorialRadius; - var e = Wgs84Eccentricity; - var s = Math.Sqrt(Math.Pow(1 + e, 1 + e) * Math.Pow(1 - e, 1 - e)); - var t = Math.Tan(Math.PI / 4d - lat / 2d) / ConformalFactor(lat); - var rho = 2d * a * ScaleFactor * t / s; - - return new Point(rho * Math.Sin(lon) + FalseEasting, rho * Math.Cos(lon) + FalseNorthing); - } - - public override Location MapToLocation(Point point) - { - point.X -= FalseEasting; - point.Y -= FalseNorthing; - - var lon = Math.Atan2(point.X, point.Y); - var rho = Math.Sqrt(point.X * point.X + point.Y * point.Y); - var a = Wgs84EquatorialRadius; - var e = Wgs84Eccentricity; - var s = Math.Sqrt(Math.Pow(1 + e, 1 + e) * Math.Pow(1 - e, 1 - e)); - var t = rho * s / (2d * a * ScaleFactor); - var lat = Math.PI / 2d - 2d * Math.Atan(t); - var relChange = 1d; - - for (int i = 0; i < MaxIterations && relChange > ConvergenceTolerance; i++) - { - var newLat = Math.PI / 2d - 2d * Math.Atan(t * ConformalFactor(lat)); - relChange = Math.Abs(1d - newLat / lat); - lat = newLat; - } - - if (IsNorth) - { - lon = Math.PI - lon; - } - else - { - lat = -lat; - } - - return new Location(lat * 180d / Math.PI, lon * 180d / Math.PI); - } - - private static double ConformalFactor(double lat) - { - var eSinLat = Wgs84Eccentricity * Math.Sin(lat); - - return Math.Pow((1d - eSinLat) / (1d + eSinLat), Wgs84Eccentricity / 2d); - } - } - - /// - /// Universal Polar Stereographic North Projection - EPSG:32661. - /// - public class UpsNorthProjection : PolarStereographicProjection - { - public const string DefaultCrsId = "EPSG:32661"; - - public UpsNorthProjection() - : base(true, 0.994, 2e6, 2e6) - { - CrsId = DefaultCrsId; - } - } - - /// - /// Universal Polar Stereographic South Projection - EPSG:32761. - /// - public class UpsSouthProjection : PolarStereographicProjection - { - public const string DefaultCrsId = "EPSG:32761"; - - public UpsSouthProjection() - : base(false, 0.994, 2e6, 2e6) - { - CrsId = DefaultCrsId; - } - } -} diff --git a/MapProjections/Shared/WebMercatorProjection.cs b/MapProjections/Shared/WebMercatorProjection.cs index b10bb785..19b4f559 100644 --- a/MapProjections/Shared/WebMercatorProjection.cs +++ b/MapProjections/Shared/WebMercatorProjection.cs @@ -4,9 +4,6 @@ using ProjNet.CoordinateSystems; using System; -#if !UWP -using System.Windows; -#endif namespace MapControl.Projections { diff --git a/MapProjections/Shared/AutoUtmProjection.cs b/MapProjections/Shared/Wgs84AutoUtmProjection.cs similarity index 53% rename from MapProjections/Shared/AutoUtmProjection.cs rename to MapProjections/Shared/Wgs84AutoUtmProjection.cs index 8fd1d0a8..3b7c00cf 100644 --- a/MapProjections/Shared/AutoUtmProjection.cs +++ b/MapProjections/Shared/Wgs84AutoUtmProjection.cs @@ -4,50 +4,47 @@ using ProjNet.CoordinateSystems; using System; -using System.Windows; namespace MapControl.Projections { - public class AutoUtmProjection : GeoApiProjection + public class Wgs84AutoUtmProjection : GeoApiProjection { public const string DefaultCrsId = "AUTO2:42001"; - public int ZoneNumber { get; private set; } - public bool ZoneIsNorth { get; private set; } - - public AutoUtmProjection() + public Wgs84AutoUtmProjection(bool useZoneCrsId = false) { + UseZoneCrsId = useZoneCrsId; UpdateZone(); } - public bool UseZoneCrsId { get; set; } + public int Zone { get; private set; } + public bool IsNorth { get; private set; } + public bool UseZoneCrsId { get; } - public override Point? LocationToMap(Location location) + public override Location Center { - UpdateZone(); - - return base.LocationToMap(location); - } - - public override Location MapToLocation(Point point) - { - UpdateZone(); - - return base.MapToLocation(point); + get => base.Center; + set + { + if (!Equals(base.Center, value)) + { + base.Center = value; + UpdateZone(); + } + } } private void UpdateZone() { - var north = Center.Latitude >= 0d; var lon = Location.NormalizeLongitude(Center.Longitude); var zone = (int)Math.Floor(lon / 6d) + 31; + var north = Center.Latitude >= 0d; - if (ZoneNumber != zone || ZoneIsNorth != north) + if (Zone != zone || IsNorth != north) { - ZoneNumber = zone; - ZoneIsNorth = north; - - CoordinateSystem = ProjectedCoordinateSystem.WGS84_UTM(ZoneNumber, ZoneIsNorth); + Zone = zone; + IsNorth = north; + CoordinateSystem = ProjectedCoordinateSystem.WGS84_UTM(Zone, IsNorth); if (!UseZoneCrsId) { diff --git a/MapProjections/Shared/WorldMercatorProjection.cs b/MapProjections/Shared/WorldMercatorProjection.cs index 3a6a158e..185f081e 100644 --- a/MapProjections/Shared/WorldMercatorProjection.cs +++ b/MapProjections/Shared/WorldMercatorProjection.cs @@ -3,9 +3,6 @@ // Licensed under the Microsoft Public License (Ms-PL) using System; -#if !UWP -using System.Windows; -#endif namespace MapControl.Projections { diff --git a/MapProjections/UWP/MapProjections.UWP.csproj b/MapProjections/UWP/MapProjections.UWP.csproj index 864540c7..d09d81ed 100644 --- a/MapProjections/UWP/MapProjections.UWP.csproj +++ b/MapProjections/UWP/MapProjections.UWP.csproj @@ -40,9 +40,6 @@ PackageReference - - AutoUtmProjection.cs - Ed50UtmProjection.cs @@ -55,12 +52,12 @@ GeoApiProjectionFactory.cs - - PolarStereographicProjection.cs - WebMercatorProjection.cs + + Wgs84AutoUtmProjection.cs + Wgs84UtmProjection.cs