ClemensFischer/XAML-Map-Control
1
0
Fork 0
mirror of https://github.com/ClemensFischer/XAML-Map-Control.git synced 2026-01-08 01:30:56 +01:00
Code Issues Actions Packages Projects Releases Wiki Activity

Add map projections

Browse source 
Add TransverseMercatorProjection, PolarStereographicProjection and derived UTM/UPS projections to the base library.
This commit is contained in:
ClemensFischer 2022-12-13 18:22:18 +01:00
parent 568f55feb0
commit e45180b26a
18 changed files with 647 additions and 332 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
MapControl/Shared/EquirectangularProjection.cs
View file

@ -17,7 +17,8 @@ namespace MapControl
/// </summary>
public class EquirectangularProjection : MapProjection
{
public const string DefaultCrsId = "EPSG:4326";
public const int DefaultEpsgCode = 4326;
public static readonly string DefaultCrsId = $"EPSG:{DefaultEpsgCode}";
public EquirectangularProjection()
{

13
MapControl/Shared/MapProjection.cs
View file

@ -44,15 +44,12 @@ namespace MapControl
/// <summary>
/// Gets or sets an optional projection center.
/// </summary>
public Location Center { get; set; } = new Location();
public virtual Location Center { get; set; } = new Location();
/// <summary>
/// Gets the relative map scale at the specified Location.
/// </summary>
public virtual Scale GetRelativeScale(Location location)
{
return new Scale(1d, 1d);
}
public virtual Scale GetRelativeScale(Location location) => new Scale(1d, 1d);
/// <summary>
/// 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;
}
}
}

64
MapControl/Shared/MapProjectionFactory.cs
View file

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

144
MapControl/Shared/PolarStereographicProjection.cs Normal file
View file

@ -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
{
/// <summary>
/// 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.
/// </summary>
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);
}
}
/// <summary>
/// Universal Polar Stereographic North Projection - EPSG:32661.
/// </summary>
public class UpsNorthProjection : PolarStereographicProjection
{
public const int DefaultEpsgCode = 32661;
public static readonly string DefaultCrsId = $"EPSG:{DefaultEpsgCode}";
public UpsNorthProjection()
{
CrsId = DefaultCrsId;
IsNorth = true;
}
}
/// <summary>
/// Universal Polar Stereographic South Projection - EPSG:32761.
/// </summary>
public class UpsSouthProjection : PolarStereographicProjection
{
public const int DefaultEpsgCode = 32761;
public static readonly string DefaultCrsId = $"EPSG:{DefaultEpsgCode}";
public UpsSouthProjection()
{
CrsId = DefaultCrsId;
IsNorth = false;
}
}
}

169
MapControl/Shared/TransverseMercatorProjection.cs Normal file
View file

@ -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
{
/// <summary>
/// Universal Transverse Mercator Projection.
/// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/pp/1395/report.pdf), p.57-64.
/// </summary>
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);
}
}
}

132
MapControl/Shared/UtmProjection.cs Normal file
View file

@ -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
{
/// <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;
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;
}
}
}
/// <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;
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;
}
}
/// <summary>
/// WGS84 UTM Projection with automatic zone selection from projection center.
/// </summary>
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);
}
}
}
}

3
MapControl/Shared/WebMercatorProjection.cs
View file

@ -15,7 +15,8 @@ namespace MapControl
/// </summary>
public class WebMercatorProjection : MapProjection
{
public const string DefaultCrsId = "EPSG:3857";
public const int DefaultEpsgCode = 3857;
public static readonly string DefaultCrsId = $"EPSG:{DefaultEpsgCode}";
public WebMercatorProjection()
{

36
MapControl/Shared/WorldMercatorProjection.cs
View file

@ -15,10 +15,8 @@ namespace MapControl
/// </summary>
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)
}
}
}

9
MapControl/UWP/MapControl.UWP.csproj
View file

@ -146,6 +146,9 @@
<Compile Include="..\Shared\OrthographicProjection.cs">
<Link>OrthographicProjection.cs</Link>
</Compile>
<Compile Include="..\Shared\PolarStereographicProjection.cs">
<Link>PolarStereographicProjection.cs</Link>
</Compile>
<Compile Include="..\Shared\PushpinBorder.cs">
<Link>PushpinBorder.cs</Link>
</Compile>
@ -170,6 +173,12 @@
<Compile Include="..\Shared\TileSource.cs">
<Link>TileSource.cs</Link>
</Compile>
<Compile Include="..\Shared\TransverseMercatorProjection.cs">
<Link>TransverseMercatorProjection.cs</Link>
</Compile>
<Compile Include="..\Shared\UtmProjection.cs">
<Link>UtmProjection.cs</Link>
</Compile>
<Compile Include="..\Shared\ViewportChangedEventArgs.cs">
<Link>ViewportChangedEventArgs.cs</Link>
</Compile>

10
MapProjections/Shared/Ed50UtmProjection.cs
View file

@ -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}\"]]";
}
}
}

10
MapProjections/Shared/Etrs89UtmProjection.cs
View file

@ -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}\"]]";
}
}
}

60
MapProjections/Shared/GeoApiProjection.cs
View file

@ -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
/// </summary>
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;
}
/// <summary>
/// Gets or sets the ICoordinateSystem of the MapProjection.
/// </summary>
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);
}
}
}

124
MapProjections/Shared/GeoApiProjectionFactory.cs
View file

@ -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<int, string> CoordinateSystemWkts { get; } = new Dictionary<int, string>();
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);

142
MapProjections/Shared/PolarStereographicProjection.cs
View file

@ -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
{
/// <summary>
/// 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.
/// </summary>
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);
}
}
/// <summary>
/// Universal Polar Stereographic North Projection - EPSG:32661.
/// </summary>
public class UpsNorthProjection : PolarStereographicProjection
{
public const string DefaultCrsId = "EPSG:32661";
public UpsNorthProjection()
: base(true, 0.994, 2e6, 2e6)
{
CrsId = DefaultCrsId;
}
}
/// <summary>
/// Universal Polar Stereographic South Projection - EPSG:32761.
/// </summary>
public class UpsSouthProjection : PolarStereographicProjection
{
public const string DefaultCrsId = "EPSG:32761";
public UpsSouthProjection()
: base(false, 0.994, 2e6, 2e6)
{
CrsId = DefaultCrsId;
}
}
}

3
MapProjections/Shared/WebMercatorProjection.cs
View file

@ -4,9 +4,6 @@
using ProjNet.CoordinateSystems;
using System;
#if !UWP
using System.Windows;
#endif
namespace MapControl.Projections
{

45
MapProjections/Shared/AutoUtmProjection.cs → MapProjections/Shared/Wgs84AutoUtmProjection.cs
View file

@ -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)
{

3
MapProjections/Shared/WorldMercatorProjection.cs
View file

@ -3,9 +3,6 @@
// Licensed under the Microsoft Public License (Ms-PL)
using System;
#if !UWP
using System.Windows;
#endif
namespace MapControl.Projections
{

9
MapProjections/UWP/MapProjections.UWP.csproj
View file

@ -40,9 +40,6 @@
<RestoreProjectStyle>PackageReference</RestoreProjectStyle>
</PropertyGroup>
<ItemGroup>
<Compile Include="..\Shared\AutoUtmProjection.cs">
<Link>AutoUtmProjection.cs</Link>
</Compile>
<Compile Include="..\Shared\Ed50UtmProjection.cs">
<Link>Ed50UtmProjection.cs</Link>
</Compile>
@ -55,12 +52,12 @@
<Compile Include="..\Shared\GeoApiProjectionFactory.cs">
<Link>GeoApiProjectionFactory.cs</Link>
</Compile>
<Compile Include="..\Shared\PolarStereographicProjection.cs">
<Link>PolarStereographicProjection.cs</Link>
</Compile>
<Compile Include="..\Shared\WebMercatorProjection.cs">
<Link>WebMercatorProjection.cs</Link>
</Compile>
<Compile Include="..\Shared\Wgs84AutoUtmProjection.cs">
<Link>Wgs84AutoUtmProjection.cs</Link>
</Compile>
<Compile Include="..\Shared\Wgs84UtmProjection.cs">
<Link>Wgs84UtmProjection.cs</Link>
</Compile>