diff --git a/MapControl/Shared/Etrs89UtmProjection.cs b/MapControl/Shared/Etrs89UtmProjection.cs index e77f3809..17490e15 100644 --- a/MapControl/Shared/Etrs89UtmProjection.cs +++ b/MapControl/Shared/Etrs89UtmProjection.cs @@ -14,7 +14,7 @@ namespace MapControl public int Zone { get; } - public Etrs89UtmProjection(int zone) + public Etrs89UtmProjection(int zone) : base(zone) { if (zone < FirstZone || zone > LastZone) { @@ -27,7 +27,6 @@ namespace MapControl // GRS 1980 EquatorialRadius = 6378137d; Flattening = 1d / 298.257222101; - CentralMeridian = zone * 6d - 183d; } } } diff --git a/MapControl/Shared/MapProjection.cs b/MapControl/Shared/MapProjection.cs index f15c8d4d..29687021 100644 --- a/MapControl/Shared/MapProjection.cs +++ b/MapControl/Shared/MapProjection.cs @@ -43,16 +43,15 @@ namespace MapControl /// public string CrsId { get; protected set; } - /// - /// Indicates whether the projection is normal cylindrical, see - /// https://en.wikipedia.org/wiki/Map_projection#Normal_cylindrical. - /// + public double EquatorialRadius { get; protected set; } = Wgs84EquatorialRadius; + public double Flattening { get; protected set; } = Wgs84Flattening; + public double ScaleFactor { get; protected set; } = 1d; + public double CentralMeridian { get; protected set; } + public double LatitudeOfOrigin { get; protected set; } + public double FalseEasting { get; protected set; } + public double FalseNorthing { get; protected set; } public bool IsNormalCylindrical { get; protected set; } - public double EquatorialRadius { get; set; } = Wgs84EquatorialRadius; - public double CentralMeridian { get; set; } - public double LatitudeOfOrigin { get; set; } - /// /// Gets the grid convergence angle in degrees at the specified geographic coordinates. /// Used for rotating the Rect resulting from BoundingBoxToMap in non-normal-cylindrical @@ -64,7 +63,12 @@ namespace MapControl /// Gets the relative transform at the specified geographic coordinates. /// The returned Matrix represents the local relative scale and rotation. /// - public abstract Matrix RelativeTransform(double latitude, double longitude); + public virtual Matrix RelativeTransform(double latitude, double longitude) + { + var transform = new Matrix(ScaleFactor, 0d, 0d, ScaleFactor, 0d, 0d); + transform.Rotate(-GridConvergence(latitude, longitude)); + return transform; + } /// /// Transforms geographic coordinates to a Point in projected map coordinates. diff --git a/MapControl/Shared/MapProjectionFactory.cs b/MapControl/Shared/MapProjectionFactory.cs index 5a2a1c40..6ca6dcf1 100644 --- a/MapControl/Shared/MapProjectionFactory.cs +++ b/MapControl/Shared/MapProjectionFactory.cs @@ -59,9 +59,9 @@ namespace MapControl var c when c is >= Nad83UtmProjection.FirstZoneEpsgCode and <= Nad83UtmProjection.LastZoneEpsgCode => new Nad83UtmProjection(c % 100), var c when c is >= Wgs84UtmProjection.FirstZoneNorthEpsgCode - and <= Wgs84UtmProjection.LastZoneNorthEpsgCode => new Wgs84UtmProjection(c % 100, Hemisphere.North), + and <= Wgs84UtmProjection.LastZoneNorthEpsgCode => new Wgs84UtmProjection(c % 100, true), var c when c is >= Wgs84UtmProjection.FirstZoneSouthEpsgCode - and <= Wgs84UtmProjection.LastZoneSouthEpsgCode => new Wgs84UtmProjection(c % 100, Hemisphere.South), + and <= Wgs84UtmProjection.LastZoneSouthEpsgCode => new Wgs84UtmProjection(c % 100, false), _ => null }; } diff --git a/MapControl/Shared/Nad27UtmProjection.cs b/MapControl/Shared/Nad27UtmProjection.cs index 3cac3657..2fd35922 100644 --- a/MapControl/Shared/Nad27UtmProjection.cs +++ b/MapControl/Shared/Nad27UtmProjection.cs @@ -14,7 +14,7 @@ namespace MapControl public int Zone { get; } - public Nad27UtmProjection(int zone) + public Nad27UtmProjection(int zone) : base(zone) { if (zone < FirstZone || zone > LastZone) { @@ -27,7 +27,6 @@ namespace MapControl // Clarke 1866 EquatorialRadius = 6378206.4; Flattening = 1d / 294.978698213898; - CentralMeridian = zone * 6d - 183d; } } } diff --git a/MapControl/Shared/Nad83UtmProjection.cs b/MapControl/Shared/Nad83UtmProjection.cs index 83b3eb8d..e46566da 100644 --- a/MapControl/Shared/Nad83UtmProjection.cs +++ b/MapControl/Shared/Nad83UtmProjection.cs @@ -14,7 +14,7 @@ namespace MapControl public int Zone { get; } - public Nad83UtmProjection(int zone) + public Nad83UtmProjection(int zone) : base(zone) { if (zone < FirstZone || zone > LastZone) { @@ -27,7 +27,6 @@ namespace MapControl // GRS 1980 EquatorialRadius = 6378137d; Flattening = 1d / 298.257222101; - CentralMeridian = zone * 6d - 183d; } } } diff --git a/MapControl/Shared/PolarStereographicProjection.cs b/MapControl/Shared/PolarStereographicProjection.cs index 320ce833..c3a9ee34 100644 --- a/MapControl/Shared/PolarStereographicProjection.cs +++ b/MapControl/Shared/PolarStereographicProjection.cs @@ -15,10 +15,13 @@ namespace MapControl /// public class PolarStereographicProjection : MapProjection { - 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 PolarStereographicProjection(bool north) + { + LatitudeOfOrigin = north ? 90d : -90d; + ScaleFactor = 0.994; + FalseEasting = 2e6; + FalseNorthing = 2e6; + } public override double GridConvergence(double latitude, double longitude) { @@ -29,20 +32,17 @@ namespace MapControl { var sign = Math.Sign(LatitudeOfOrigin); var phi = sign * latitude * Math.PI / 180d; - var e = Math.Sqrt((2d - Flattening) * Flattening); var eSinPhi = e * Math.Sin(phi); - var t = Math.Tan(Math.PI / 4d - phi / 2d) / Math.Pow((1d - eSinPhi) / (1d + eSinPhi), e / 2d); // p.161 (15-9) - // r == ρ/a var r = 2d * ScaleFactor * t / Math.Sqrt(Math.Pow(1d + e, 1d + e) * Math.Pow(1d - e, 1d - e)); // p.161 (21-33) var m = Math.Cos(phi) / Math.Sqrt(1d - eSinPhi * eSinPhi); // p.160 (14-15) var k = r / m; // p.161 (21-32) var transform = new Matrix(k, 0d, 0d, k, 0d, 0d); - transform.Rotate(-sign * longitude); + transform.Rotate(-sign * (longitude - CentralMeridian)); return transform; } @@ -99,10 +99,9 @@ namespace MapControl { } - public Wgs84UpsNorthProjection(string crsId) + public Wgs84UpsNorthProjection(string crsId) : base(true) { CrsId = crsId; - LatitudeOfOrigin = 90d; } } @@ -118,10 +117,9 @@ namespace MapControl { } - public Wgs84UpsSouthProjection(string crsId) + public Wgs84UpsSouthProjection(string crsId) : base(false) { CrsId = crsId; - LatitudeOfOrigin = -90d; } } } diff --git a/MapControl/Shared/StereographicProjection.cs b/MapControl/Shared/StereographicProjection.cs index 3800c6d6..7bcae21a 100644 --- a/MapControl/Shared/StereographicProjection.cs +++ b/MapControl/Shared/StereographicProjection.cs @@ -25,7 +25,7 @@ namespace MapControl "{0},1,{1:0.########},{2:0.########}", crsId, centerLongitude, centerLatitude); } - public override double GridConvergence(double latitude, double longitude) + private void GetScaleAndGridConvergence(double latitude, double longitude, out double scale, out double gamma) { var phi0 = LatitudeOfOrigin * Math.PI / 180d; // φ1 var phi1 = latitude * Math.PI / 180d; @@ -44,13 +44,26 @@ namespace MapControl var dCosPhi = k2 * cosPhi2 - k1 * cosPhi1; var dSinPhi = k2 * sinPhi2 - k1 * sinPhi1; - return Math.Atan2(-sinLambda * dCosPhi, + scale = k1; + gamma = Math.Atan2(-sinLambda * dCosPhi, cosPhi0 * dSinPhi - sinPhi0 * cosLambda * dCosPhi) * 180d / Math.PI; } + public override double GridConvergence(double latitude, double longitude) + { + GetScaleAndGridConvergence(latitude, longitude, out double _, out double gamma); + + return gamma; + } + public override Matrix RelativeTransform(double latitude, double longitude) { - return new Matrix(1d, 0d, 0d, 1d, 0d, 0d); + GetScaleAndGridConvergence(latitude, longitude, out double scale, out double gamma); + + var transform = new Matrix(scale, 0d, 0d, scale, 0d, 0d); + transform.Rotate(-gamma); + + return transform; } public override Point LocationToMap(double latitude, double longitude) diff --git a/MapControl/Shared/TransverseMercatorProjection.cs b/MapControl/Shared/TransverseMercatorProjection.cs index 0a7c13eb..e0c0430a 100644 --- a/MapControl/Shared/TransverseMercatorProjection.cs +++ b/MapControl/Shared/TransverseMercatorProjection.cs @@ -16,6 +16,8 @@ namespace MapControl /// public class TransverseMercatorProjection : MapProjection { + private double f; + private double n; private double a1; // α1 private double a2; // α2 private double a3; // α3 @@ -28,38 +30,36 @@ namespace MapControl private double f1; // A/a private double f2; // 2*sqrt(n)/(1+n) + private void InitializeParameters() + { + f = Flattening; + n = f / (2d - f); + var n2 = n * n; + var n3 = n * n2; + a1 = n / 2d - n2 * 2d / 3d + n3 * 5d / 16d; + a2 = n2 * 13d / 48d - n3 * 3d / 5d; + a3 = n3 * 61d / 240d; + b1 = n / 2d - n2 * 2d / 3d + n3 * 37d / 96d; + b2 = n2 / 48d + n3 / 15d; + b3 = n3 * 17d / 480d; + d1 = n * 2d - n2 * 2d / 3d - n3 * 2d; + d2 = n2 * 7d / 3d - n3 * 8d / 5d; + d3 = n3 * 56d / 15d; + f1 = (1d + n2 / 4d + n2 * n2 / 64d) / (1d + n); + f2 = 2d * Math.Sqrt(n) / (1d + n); + } + public TransverseMercatorProjection() { - Flattening = Wgs84Flattening; } - public double Flattening + public TransverseMercatorProjection(int utmZone) : this() { - get; - set - { - field = value; - var n = field / (2d - field); - var n2 = n * n; - var n3 = n * n2; - a1 = n / 2d - n2 * 2d / 3d + n3 * 5d / 16d; - a2 = n2 * 13d / 48d - n3 * 3d / 5d; - a3 = n3 * 61d / 240d; - b1 = n / 2d - n2 * 2d / 3d + n3 * 37d / 96d; - b2 = n2 / 48d + n3 / 15d; - b3 = n3 * 17d / 480d; - d1 = n * 2d - n2 * 2d / 3d - n3 * 2d; - d2 = n2 * 7d / 3d - n3 * 8d / 5d; - d3 = n3 * 56d / 15d; - f1 = (1d + n2 / 4d + n2 * n2 / 64d) / (1d + n); - f2 = 2d * Math.Sqrt(n) / (1d + n); - } + CentralMeridian = utmZone * 6d - 183d; + ScaleFactor = 0.9996; + FalseEasting = 5e5; } - public double ScaleFactor { get; set; } = 0.9996; - public double FalseEasting { get; set; } = 5e5; - public double FalseNorthing { get; set; } - public override double GridConvergence(double latitude, double longitude) { // φ @@ -74,6 +74,11 @@ namespace MapControl public override Matrix RelativeTransform(double latitude, double longitude) { + if (f != Flattening) + { + InitializeParameters(); + } + // φ var phi = latitude * Math.PI / 180d; var sinPhi = Math.Sin(phi); @@ -99,7 +104,6 @@ namespace MapControl 4d * a2 * Math.Sin(4d * xi_) * Math.Sinh(4d * eta_) + 6d * a3 * Math.Sin(6d * xi_) * Math.Sinh(6d * eta_); - var n = Flattening / (2d - Flattening); var m = (1d - n) / (1d + n) * Math.Tan(phi); var k = ScaleFactor * f1 * Math.Sqrt((1d + m * m) * (sigma * sigma + tau * tau) / (t * t + cosLambda * cosLambda)); @@ -114,6 +118,11 @@ namespace MapControl public override Point LocationToMap(double latitude, double longitude) { + if (f != Flattening) + { + InitializeParameters(); + } + // φ var phi = latitude * Math.PI / 180d; var sinPhi = Math.Sin(phi); @@ -143,6 +152,11 @@ namespace MapControl public override Location MapToLocation(double x, double y) { + if (f != Flattening) + { + InitializeParameters(); + } + // k0 * A var k0A = ScaleFactor * EquatorialRadius * f1; // ξ diff --git a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs index de62674e..12026aeb 100644 --- a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs +++ b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs @@ -15,10 +15,16 @@ namespace MapControl /// public class TransverseMercatorProjectionSnyder : MapProjection { - public double Flattening { get; set; } = Wgs84Flattening; - public double ScaleFactor { get; set; } = 0.9996; - public double FalseEasting { get; set; } = 5e5; - public double FalseNorthing { get; set; } + public TransverseMercatorProjectionSnyder() + { + } + + public TransverseMercatorProjectionSnyder(int utmZone) : this() + { + CentralMeridian = utmZone * 6d - 183d; + ScaleFactor = 0.9996; + FalseEasting = 5e5; + } public override double GridConvergence(double latitude, double longitude) { diff --git a/MapControl/Shared/Wgs84UtmProjection.cs b/MapControl/Shared/Wgs84UtmProjection.cs index 470d859f..82aae18a 100644 --- a/MapControl/Shared/Wgs84UtmProjection.cs +++ b/MapControl/Shared/Wgs84UtmProjection.cs @@ -2,12 +2,6 @@ namespace MapControl { - public enum Hemisphere - { - North, - South - } - /// /// WGS84 Universal Transverse Mercator Projection - /// EPSG:32601 to EPSG:32660 and EPSG:32701 to EPSG:32760. @@ -22,9 +16,8 @@ namespace MapControl public const int LastZoneSouthEpsgCode = 32700 + LastZone; public int Zone { get; } - public Hemisphere Hemisphere { get; } - public Wgs84UtmProjection(int zone, Hemisphere hemisphere) + public Wgs84UtmProjection(int zone, bool north) : base(zone) { if (zone < FirstZone || zone > LastZone) { @@ -32,10 +25,8 @@ namespace MapControl } Zone = zone; - Hemisphere = hemisphere; - CrsId = $"EPSG:{(hemisphere == Hemisphere.North ? 32600 : 32700) + zone}"; - CentralMeridian = zone * 6d - 183d; - FalseNorthing = hemisphere == Hemisphere.North ? 0d : 1e7; + CrsId = $"EPSG:{(north ? 32600 : 32700) + zone}"; + FalseNorthing = north ? 0d : 1e7; } } } diff --git a/MapControl/Shared/WorldMercatorProjection.cs b/MapControl/Shared/WorldMercatorProjection.cs index 631f5bb7..c1560320 100644 --- a/MapControl/Shared/WorldMercatorProjection.cs +++ b/MapControl/Shared/WorldMercatorProjection.cs @@ -16,8 +16,6 @@ namespace MapControl { public const string DefaultCrsId = "EPSG:3395"; - public double Flattening { get; set; } = Wgs84Flattening; - public WorldMercatorProjection() // parameterless constructor for XAML : this(DefaultCrsId) { @@ -31,8 +29,8 @@ namespace MapControl public override Matrix RelativeTransform(double latitude, double longitude) { - var phi = latitude * Math.PI / 180d; var e2 = (2d - Flattening) * Flattening; + var phi = latitude * Math.PI / 180d; var sinPhi = Math.Sin(phi); var k = Math.Sqrt(1d - e2 * sinPhi * sinPhi) / Math.Cos(phi); // p.44 (7-8) diff --git a/MapProjections/Shared/Ed50UtmProjection.cs b/MapProjections/Shared/Ed50UtmProjection.cs deleted file mode 100644 index 0ef54a3f..00000000 --- a/MapProjections/Shared/Ed50UtmProjection.cs +++ /dev/null @@ -1,46 +0,0 @@ -using System; - -namespace MapControl.Projections -{ - /// - /// ED50 Universal Transverse Mercator Projection. - /// - public class Ed50UtmProjection : ProjNetMapProjection - { - public const int FirstZone = 28; - public const int LastZone = 38; - public const int FirstZoneEpsgCode = 23000 + FirstZone; - public const int LastZoneEpsgCode = 23000 + LastZone; - - public int Zone { get; } - - public Ed50UtmProjection(int zone) - : base(new TransverseMercatorProjection - { - EquatorialRadius = 6378388d, - Flattening = 297, - CentralMeridian = 6 * zone - 183, - }) - { - if (zone < FirstZone || zone > LastZone) - { - throw new ArgumentException($"Invalid ED50 UTM zone {zone}.", nameof(zone)); - } - - Zone = zone; - CoordinateSystemWkt = - $"PROJCS[\"ED50 / UTM zone {zone}N\"," + - WktConstants.GeogCsEd50 + "," + - "PROJECTION[\"Transverse_Mercator\"]," + - "PARAMETER[\"latitude_of_origin\",0]," + - $"PARAMETER[\"central_meridian\",{6 * zone - 183}]," + - "PARAMETER[\"scale_factor\",0.9996]," + - "PARAMETER[\"false_easting\",500000]," + - "PARAMETER[\"false_northing\",0]," + - "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + - "AXIS[\"Easting\",EAST]," + - "AXIS[\"Northing\",NORTH]," + - $"AUTHORITY[\"EPSG\",\"230{zone:00}\"]]"; - } - } -} diff --git a/MapProjections/Shared/Etrs89UtmProjection.cs b/MapProjections/Shared/Etrs89UtmProjection.cs deleted file mode 100644 index e1b21514..00000000 --- a/MapProjections/Shared/Etrs89UtmProjection.cs +++ /dev/null @@ -1,29 +0,0 @@ -namespace MapControl.Projections -{ - /// - /// ETRS89 Universal Transverse Mercator Projection - EPSG:25828 to EPSG:25838. - /// - public class Etrs89UtmProjection : ProjNetMapProjection - { - public int Zone { get; } - - public Etrs89UtmProjection(int zone) - : base(new MapControl.Etrs89UtmProjection(zone)) - { - Zone = zone; - CoordinateSystemWkt = - $"PROJCS[\"ETRS89 / UTM zone {zone}N\"," + - WktConstants.GeogCsEtrs89 + "," + - "PROJECTION[\"Transverse_Mercator\"]," + - "PARAMETER[\"latitude_of_origin\",0]," + - $"PARAMETER[\"central_meridian\",{6 * zone - 183}]," + - "PARAMETER[\"scale_factor\",0.9996]," + - "PARAMETER[\"false_easting\",500000]," + - "PARAMETER[\"false_northing\",0]," + - "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + - "AXIS[\"Easting\",EAST]," + - "AXIS[\"Northing\",NORTH]," + - $"AUTHORITY[\"EPSG\",\"258{zone:00}\"]]"; - } - } -} diff --git a/MapProjections/Shared/Nad27UtmProjection.cs b/MapProjections/Shared/Nad27UtmProjection.cs deleted file mode 100644 index 098b5e42..00000000 --- a/MapProjections/Shared/Nad27UtmProjection.cs +++ /dev/null @@ -1,29 +0,0 @@ -namespace MapControl.Projections -{ - /// - /// NAD27 Universal Transverse Mercator Projection - EPSG:26701 to EPSG:26722. - /// - public class Nad27UtmProjection : ProjNetMapProjection - { - public int Zone { get; } - - public Nad27UtmProjection(int zone) - : base(new MapControl.Nad27UtmProjection(zone)) - { - Zone = zone; - CoordinateSystemWkt = - $"PROJCS[\"NAD27 / UTM zone {zone}N\"," + - WktConstants.GeogCsNad27 + "," + - "PROJECTION[\"Transverse_Mercator\"]," + - "PARAMETER[\"latitude_of_origin\",0]," + - $"PARAMETER[\"central_meridian\",{6 * zone - 183}]," + - "PARAMETER[\"scale_factor\",0.9996]," + - "PARAMETER[\"false_easting\",500000]," + - "PARAMETER[\"false_northing\",0]," + - "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + - "AXIS[\"Easting\",EAST]," + - "AXIS[\"Northing\",NORTH]," + - $"AUTHORITY[\"EPSG\",\"267{zone:00}\"]]"; - } - } -} diff --git a/MapProjections/Shared/Nad83UtmProjection.cs b/MapProjections/Shared/Nad83UtmProjection.cs deleted file mode 100644 index 9c39dccf..00000000 --- a/MapProjections/Shared/Nad83UtmProjection.cs +++ /dev/null @@ -1,29 +0,0 @@ -namespace MapControl.Projections -{ - /// - /// NAD83 Universal Transverse Mercator Projection - EPSG:26901 to EPSG:26923. - /// - public class Nad83UtmProjection : ProjNetMapProjection - { - public int Zone { get; } - - public Nad83UtmProjection(int zone) - : base(new MapControl.Nad83UtmProjection(zone)) - { - Zone = zone; - CoordinateSystemWkt = - $"PROJCS[\"NAD83 / UTM zone {zone}N\"," + - WktConstants.GeogCsNad83 + "," + - "PROJECTION[\"Transverse_Mercator\"]," + - "PARAMETER[\"latitude_of_origin\",0]," + - $"PARAMETER[\"central_meridian\",{6 * zone - 183}]," + - "PARAMETER[\"scale_factor\",0.9996]," + - "PARAMETER[\"false_easting\",500000]," + - "PARAMETER[\"false_northing\",0]," + - "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + - "AXIS[\"Easting\",EAST]," + - "AXIS[\"Northing\",NORTH]," + - $"AUTHORITY[\"EPSG\",\"269{zone:00}\"]]"; - } - } -} diff --git a/MapProjections/Shared/ProjNetMapProjection.cs b/MapProjections/Shared/ProjNetMapProjection.cs index 314bf5ba..7717668f 100644 --- a/MapProjections/Shared/ProjNetMapProjection.cs +++ b/MapProjections/Shared/ProjNetMapProjection.cs @@ -15,11 +15,9 @@ namespace MapControl.Projections /// public class ProjNetMapProjection : MapProjection { - protected MapProjection FallbackProjection { get; private set; } - - protected ProjNetMapProjection(MapProjection fallbackProjection) + public ProjNetMapProjection(ProjectedCoordinateSystem coordinateSystem) { - FallbackProjection = fallbackProjection; + CoordinateSystem = coordinateSystem; } public ProjNetMapProjection(string coordinateSystemWkt) @@ -52,7 +50,9 @@ namespace MapControl.Projections var projection = field.Projection ?? throw new ArgumentException("CoordinateSystem.Projection must not be null.", nameof(value)); - var ellipsoid = field.GeographicCoordinateSystem.HorizontalDatum.Ellipsoid; + IsNormalCylindrical = projection.Name.Contains("Pseudo-Mercator") || + projection.Name.StartsWith("Mercator"); + var transformFactory = new CoordinateTransformationFactory(); CrsId = !string.IsNullOrEmpty(field.Authority) && field.AuthorityCode > 0 @@ -67,47 +67,24 @@ namespace MapControl.Projections .CreateFromCoordinateSystems(field, GeographicCoordinateSystem.WGS84) .MathTransform; - if (projection.Name.Contains("Pseudo-Mercator")) - { - IsNormalCylindrical = true; - FallbackProjection = new MapControl.WebMercatorProjection - { - EquatorialRadius = ellipsoid.SemiMajorAxis - }; - } - else if (projection.Name.StartsWith("Mercator")) - { - IsNormalCylindrical = true; - FallbackProjection = new MapControl.WorldMercatorProjection - { - EquatorialRadius = ellipsoid.SemiMajorAxis, - Flattening = 1d / ellipsoid.InverseFlattening - }; - } - else if (projection.Name.StartsWith("Transverse_Mercator")) - { - FallbackProjection = new TransverseMercatorProjection - { - EquatorialRadius = ellipsoid.SemiMajorAxis, - Flattening = 1d / ellipsoid.InverseFlattening, - CentralMeridian = projection.GetParameter("central_meridian").Value, - ScaleFactor = projection.GetParameter("scale_factor").Value, - FalseEasting = projection.GetParameter("false_easting").Value, - FalseNorthing = projection.GetParameter("false_northing").Value - }; - } - else if (projection.Name.StartsWith("Polar_Stereographic")) - { - FallbackProjection = new PolarStereographicProjection - { - EquatorialRadius = ellipsoid.SemiMajorAxis, - Flattening = 1d / ellipsoid.InverseFlattening, - ScaleFactor = projection.GetParameter("scale_factor").Value, - FalseEasting = projection.GetParameter("false_easting").Value, - FalseNorthing = projection.GetParameter("false_northing").Value, - LatitudeOfOrigin = projection.GetParameter("latitude_of_origin").Value - }; - } + var ellipsoid = field.HorizontalDatum.Ellipsoid; + EquatorialRadius = ellipsoid.SemiMajorAxis; + Flattening = 1d / ellipsoid.InverseFlattening; + + var parameter = projection.GetParameter("scale_factor"); + ScaleFactor = parameter != null ? parameter.Value : 1d; + + parameter = projection.GetParameter("central_meridian"); + CentralMeridian = parameter != null ? parameter.Value : 0d; + + parameter = projection.GetParameter("latitude_of_origin"); + LatitudeOfOrigin = parameter != null ? parameter.Value : 0d; + + parameter = projection.GetParameter("false_easting"); + FalseEasting = parameter != null ? parameter.Value : 0d; + + parameter = projection.GetParameter("false_northing"); + FalseNorthing = parameter != null ? parameter.Value : 0d; } } @@ -139,18 +116,91 @@ namespace MapControl.Projections return new Location(coordinate[1], coordinate[0]); } - public override Matrix RelativeTransform(double latitude, double longitude) + public override double GridConvergence(double latitude, double longitude) { - return FallbackProjection != null - ? FallbackProjection.RelativeTransform(latitude, longitude) - : new Matrix(1d, 0d, 0d, 1d, 0d, 0d); + var projection = CoordinateSystem.Projection.Name; + + if (projection.StartsWith("Transverse_Mercator")) + { + return TransverseMercatorGridConvergence(latitude, longitude); + } + + if (projection.StartsWith("Polar_Stereographic")) + { + return PolarStereographicGridConvergence(longitude); + } + + return base.GridConvergence(latitude, longitude); } - public override double GridConvergence(double x, double y) + public override Matrix RelativeTransform(double latitude, double longitude) { - return FallbackProjection != null - ? FallbackProjection.GridConvergence(x, y) - : 0d; + var projection = CoordinateSystem.Projection.Name; + + if (projection.Contains("Pseudo-Mercator")) + { + return WebMercatorRelativeTransform(latitude); + } + + if (projection.StartsWith("Mercator")) + { + return WorldMercatorRelativeTransform(latitude); + } + + if (projection.StartsWith("Polar_Stereographic")) + { + return PolarStereographicRelativeTransform(latitude, longitude); + } + + return base.RelativeTransform(latitude, longitude); + } + + protected static Matrix WebMercatorRelativeTransform(double latitude) + { + var k = 1d / Math.Cos(latitude * Math.PI / 180d); // p.44 (7-3) + + return new Matrix(k, 0d, 0d, k, 0d, 0d); + } + + protected Matrix WorldMercatorRelativeTransform(double latitude) + { + var e2 = (2d - Flattening) * Flattening; + var phi = latitude * Math.PI / 180d; + var sinPhi = Math.Sin(phi); + var k = Math.Sqrt(1d - e2 * sinPhi * sinPhi) / Math.Cos(phi); // p.44 (7-8) + + return new Matrix(k, 0d, 0d, k, 0d, 0d); + } + + protected double TransverseMercatorGridConvergence(double latitude, double longitude) + { + return 180d / Math.PI * Math.Atan( + Math.Tan((longitude - CentralMeridian) * Math.PI / 180d) * + Math.Sin(latitude * Math.PI / 180d)); + } + + protected double PolarStereographicGridConvergence(double longitude) + { + return Math.Sign(LatitudeOfOrigin) * (longitude - CentralMeridian); + } + + protected Matrix PolarStereographicRelativeTransform(double latitude, double longitude) + { + var sign = Math.Sign(LatitudeOfOrigin); + var phi = sign * latitude * Math.PI / 180d; + var e = Math.Sqrt((2d - Flattening) * Flattening); + var eSinPhi = e * Math.Sin(phi); + var t = Math.Tan(Math.PI / 4d - phi / 2d) + / Math.Pow((1d - eSinPhi) / (1d + eSinPhi), e / 2d); // p.161 (15-9) + // r == ρ/a + var r = 2d * ScaleFactor * t / Math.Sqrt(Math.Pow(1d + e, 1d + e) * Math.Pow(1d - e, 1d - e)); // p.161 (21-33) + var m = Math.Cos(phi) / Math.Sqrt(1d - eSinPhi * eSinPhi); // p.160 (14-15) + var k = r / m; // p.161 (21-32) + + var transform = new Matrix(k, 0d, 0d, k, 0d, 0d); + transform.Rotate(-sign * (longitude - CentralMeridian)); + + return transform; } } } diff --git a/MapProjections/Shared/ProjNetMapProjectionFactory.cs b/MapProjections/Shared/ProjNetMapProjectionFactory.cs index aa027d6e..e1371964 100644 --- a/MapProjections/Shared/ProjNetMapProjectionFactory.cs +++ b/MapProjections/Shared/ProjNetMapProjectionFactory.cs @@ -1,4 +1,5 @@ -using System.Collections.Generic; +using ProjNet.CoordinateSystems; +using System.Collections.Generic; namespace MapControl.Projections { @@ -54,11 +55,175 @@ namespace MapControl.Projections var c when c is >= MapControl.Nad83UtmProjection.FirstZoneEpsgCode and <= MapControl.Nad83UtmProjection.LastZoneEpsgCode => new Nad83UtmProjection(c % 100), var c when c is >= MapControl.Wgs84UtmProjection.FirstZoneNorthEpsgCode - and <= MapControl.Wgs84UtmProjection.LastZoneNorthEpsgCode => new Wgs84UtmProjection(c % 100, Hemisphere.North), + and <= MapControl.Wgs84UtmProjection.LastZoneNorthEpsgCode => new Wgs84UtmProjection(c % 100, true), var c when c is >= MapControl.Wgs84UtmProjection.FirstZoneSouthEpsgCode - and <= MapControl.Wgs84UtmProjection.LastZoneSouthEpsgCode => new Wgs84UtmProjection(c % 100, Hemisphere.South), + and <= MapControl.Wgs84UtmProjection.LastZoneSouthEpsgCode => new Wgs84UtmProjection(c % 100, false), _ => base.CreateProjection(epsgCode) }; } } + /// + /// Spherical Mercator Projection - EPSG:3857, + /// implemented by setting the CoordinateSystem property of a ProjNetMapProjection. + /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/pp/1395/report.pdf), p.41-44. + /// + public class WebMercatorProjection : ProjNetMapProjection + { + public WebMercatorProjection() + : base(ProjectedCoordinateSystem.WebMercator) + { + } + } + + /// + /// Elliptical Mercator Projection - EPSG:3395, + /// implemented by setting the CoordinateSystemWkt property of a ProjNetMapProjection. + /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/pp/1395/report.pdf), p.44-45. + /// + public class WorldMercatorProjection : ProjNetMapProjection + { + public WorldMercatorProjection() : base( + "PROJCS[\"WGS 84 / World Mercator\"," + + WktConstants.GeogCsWgs84 + "," + + "PROJECTION[\"Mercator_1SP\"]," + + "PARAMETER[\"latitude_of_origin\",0]," + + "PARAMETER[\"central_meridian\",0]," + + "PARAMETER[\"scale_factor\",1]," + + "PARAMETER[\"false_easting\",0]," + + "PARAMETER[\"false_northing\",0]," + + "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + + "AXIS[\"Easting\",EAST]," + + "AXIS[\"Northing\",NORTH]," + + "AUTHORITY[\"EPSG\",\"3395\"]]") + { + } + } + + /// + /// WGS84 Universal Transverse Mercator Projection - + /// EPSG:32601 to EPSG:32660 and EPSG:32701 to EPSG:32760. + /// + public class Wgs84UtmProjection(int zone, bool north) : ProjNetMapProjection( + ProjectedCoordinateSystem.WGS84_UTM(zone, north)) + { + public int Zone => zone; + } + + /// + /// ETRS89 Universal Transverse Mercator Projection - EPSG:25828 to EPSG:25838. + /// + public class Etrs89UtmProjection(int zone) : ProjNetMapProjection( + $"PROJCS[\"ETRS89 / UTM zone {zone}N\"," + + WktConstants.GeogCsEtrs89 + "," + + "PROJECTION[\"Transverse_Mercator\"]," + + "PARAMETER[\"latitude_of_origin\",0]," + + $"PARAMETER[\"central_meridian\",{6 * zone - 183}]," + + "PARAMETER[\"scale_factor\",0.9996]," + + "PARAMETER[\"false_easting\",500000]," + + "PARAMETER[\"false_northing\",0]," + + "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + + "AXIS[\"Easting\",EAST]," + + "AXIS[\"Northing\",NORTH]," + + $"AUTHORITY[\"EPSG\",\"258{zone:00}\"]]") + { + public int Zone => zone; + } + + /// + /// NAD27 Universal Transverse Mercator Projection - EPSG:26701 to EPSG:26722. + /// + public class Nad27UtmProjection(int zone) : ProjNetMapProjection( + $"PROJCS[\"NAD27 / UTM zone {zone}N\"," + + WktConstants.GeogCsNad27 + "," + + "PROJECTION[\"Transverse_Mercator\"]," + + "PARAMETER[\"latitude_of_origin\",0]," + + $"PARAMETER[\"central_meridian\",{6 * zone - 183}]," + + "PARAMETER[\"scale_factor\",0.9996]," + + "PARAMETER[\"false_easting\",500000]," + + "PARAMETER[\"false_northing\",0]," + + "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + + "AXIS[\"Easting\",EAST]," + + "AXIS[\"Northing\",NORTH]," + + $"AUTHORITY[\"EPSG\",\"267{zone:00}\"]]") + { + public int Zone => zone; + } + /// + /// NAD83 Universal Transverse Mercator Projection - EPSG:26901 to EPSG:26923. + /// + public class Nad83UtmProjection(int zone) : ProjNetMapProjection( + $"PROJCS[\"NAD83 / UTM zone {zone}N\"," + + WktConstants.GeogCsNad83 + "," + + "PROJECTION[\"Transverse_Mercator\"]," + + "PARAMETER[\"latitude_of_origin\",0]," + + $"PARAMETER[\"central_meridian\",{6 * zone - 183}]," + + "PARAMETER[\"scale_factor\",0.9996]," + + "PARAMETER[\"false_easting\",500000]," + + "PARAMETER[\"false_northing\",0]," + + "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + + "AXIS[\"Easting\",EAST]," + + "AXIS[\"Northing\",NORTH]," + + $"AUTHORITY[\"EPSG\",\"269{zone:00}\"]]") + { + public int Zone => zone; + } + + /// + /// ED50 Universal Transverse Mercator Projection. + /// + public class Ed50UtmProjection(int zone) : ProjNetMapProjection( + $"PROJCS[\"ED50 / UTM zone {zone}N\"," + + WktConstants.GeogCsEd50 + "," + + "PROJECTION[\"Transverse_Mercator\"]," + + "PARAMETER[\"latitude_of_origin\",0]," + + $"PARAMETER[\"central_meridian\",{6 * zone - 183}]," + + "PARAMETER[\"scale_factor\",0.9996]," + + "PARAMETER[\"false_easting\",500000]," + + "PARAMETER[\"false_northing\",0]," + + "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + + "AXIS[\"Easting\",EAST]," + + "AXIS[\"Northing\",NORTH]," + + $"AUTHORITY[\"EPSG\",\"230{zone:00}\"]]") + { + public const int FirstZone = 28; + public const int LastZone = 38; + public const int FirstZoneEpsgCode = 23000 + FirstZone; + public const int LastZoneEpsgCode = 23000 + LastZone; + + public int Zone { get; } = zone; + } + + public class Wgs84UpsNorthProjection : ProjNetMapProjection + { + public Wgs84UpsNorthProjection() : base( + "PROJCS[\"WGS 84 / UPS North (N,E)\"," + + WktConstants.GeogCsWgs84 + "," + + "PROJECTION[\"Polar_Stereographic\"]," + + "PARAMETER[\"latitude_of_origin\",90]," + + "PARAMETER[\"central_meridian\",0]," + + "PARAMETER[\"scale_factor\",0.994]," + + "PARAMETER[\"false_easting\",2000000]," + + "PARAMETER[\"false_northing\",2000000]," + + "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + + "AUTHORITY[\"EPSG\",\"32661\"]]") + { + } + } + + public class Wgs84UpsSouthProjection : ProjNetMapProjection + { + public Wgs84UpsSouthProjection() : base( + "PROJCS[\"WGS 84 / UPS South (N,E)\"," + + WktConstants.GeogCsWgs84 + "," + + "PROJECTION[\"Polar_Stereographic\"]," + + "PARAMETER[\"latitude_of_origin\",-90]," + + "PARAMETER[\"central_meridian\",0]," + + "PARAMETER[\"scale_factor\",0.994]," + + "PARAMETER[\"false_easting\",2000000]," + + "PARAMETER[\"false_northing\",2000000]," + + "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + + "AUTHORITY[\"EPSG\",\"32761\"]]") + { + } + } } diff --git a/MapProjections/Shared/WebMercatorProjection.cs b/MapProjections/Shared/WebMercatorProjection.cs deleted file mode 100644 index 37659a80..00000000 --- a/MapProjections/Shared/WebMercatorProjection.cs +++ /dev/null @@ -1,17 +0,0 @@ -using ProjNet.CoordinateSystems; - -namespace MapControl.Projections -{ - /// - /// Spherical Mercator Projection implemented by setting the CoordinateSystem property of a ProjNetMapProjection. - /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/pp/1395/report.pdf), p.41-44. - /// - public class WebMercatorProjection : ProjNetMapProjection - { - public WebMercatorProjection() - : base(new MapControl.WebMercatorProjection()) - { - CoordinateSystem = ProjectedCoordinateSystem.WebMercator; - } - } -} diff --git a/MapProjections/Shared/Wgs84UpsProjections.cs b/MapProjections/Shared/Wgs84UpsProjections.cs deleted file mode 100644 index e72da80e..00000000 --- a/MapProjections/Shared/Wgs84UpsProjections.cs +++ /dev/null @@ -1,44 +0,0 @@ -#if WPF -using System.Windows.Media; -#endif - -namespace MapControl.Projections -{ - public class Wgs84UpsNorthProjection : ProjNetMapProjection - { - public Wgs84UpsNorthProjection() - : base(new MapControl.Wgs84UpsNorthProjection()) - { - CoordinateSystemWkt = - "PROJCS[\"WGS 84 / UPS North (N,E)\"," + - WktConstants.GeogCsWgs84 + "," + - "PROJECTION[\"Polar_Stereographic\"]," + - "PARAMETER[\"latitude_of_origin\",90]," + - "PARAMETER[\"central_meridian\",0]," + - "PARAMETER[\"scale_factor\",0.994]," + - "PARAMETER[\"false_easting\",2000000]," + - "PARAMETER[\"false_northing\",2000000]," + - "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + - "AUTHORITY[\"EPSG\",\"32661\"]]"; - } - } - - public class Wgs84UpsSouthProjection : ProjNetMapProjection - { - public Wgs84UpsSouthProjection() - : base(new MapControl.Wgs84UpsSouthProjection()) - { - CoordinateSystemWkt = - "PROJCS[\"WGS 84 / UPS South (N,E)\"," + - WktConstants.GeogCsWgs84 + "," + - "PROJECTION[\"Polar_Stereographic\"]," + - "PARAMETER[\"latitude_of_origin\",-90]," + - "PARAMETER[\"central_meridian\",0]," + - "PARAMETER[\"scale_factor\",0.994]," + - "PARAMETER[\"false_easting\",2000000]," + - "PARAMETER[\"false_northing\",2000000]," + - "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + - "AUTHORITY[\"EPSG\",\"32761\"]]"; - } - } -} diff --git a/MapProjections/Shared/Wgs84UtmProjection.cs b/MapProjections/Shared/Wgs84UtmProjection.cs deleted file mode 100644 index b6da1918..00000000 --- a/MapProjections/Shared/Wgs84UtmProjection.cs +++ /dev/null @@ -1,22 +0,0 @@ -using ProjNet.CoordinateSystems; - -namespace MapControl.Projections -{ - /// - /// WGS84 Universal Transverse Mercator Projection - - /// EPSG:32601 to EPSG:32660 and EPSG:32701 to EPSG:32760. - /// - public class Wgs84UtmProjection : ProjNetMapProjection - { - public int Zone { get; } - public Hemisphere Hemisphere { get; } - - public Wgs84UtmProjection(int zone, Hemisphere hemisphere) - : base(new MapControl.Wgs84UtmProjection(zone, hemisphere)) - { - Zone = zone; - Hemisphere = hemisphere; - CoordinateSystem = ProjectedCoordinateSystem.WGS84_UTM(zone, hemisphere == Hemisphere.North); - } - } -} diff --git a/MapProjections/Shared/WorldMercatorProjection.cs b/MapProjections/Shared/WorldMercatorProjection.cs deleted file mode 100644 index 77a512ca..00000000 --- a/MapProjections/Shared/WorldMercatorProjection.cs +++ /dev/null @@ -1,27 +0,0 @@ -namespace MapControl.Projections -{ - /// - /// Elliptical Mercator Projection implemented by setting the WKT property of a ProjNetMapProjection. - /// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/pp/1395/report.pdf), p.44-45. - /// - public class WorldMercatorProjection : ProjNetMapProjection - { - public WorldMercatorProjection() - : base(new MapControl.WorldMercatorProjection()) - { - CoordinateSystemWkt = - "PROJCS[\"WGS 84 / World Mercator\"," + - WktConstants.GeogCsWgs84 + "," + - "PROJECTION[\"Mercator_1SP\"]," + - "PARAMETER[\"latitude_of_origin\",0]," + - "PARAMETER[\"central_meridian\",0]," + - "PARAMETER[\"scale_factor\",1]," + - "PARAMETER[\"false_easting\",0]," + - "PARAMETER[\"false_northing\",0]," + - "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + - "AXIS[\"Easting\",EAST]," + - "AXIS[\"Northing\",NORTH]," + - "AUTHORITY[\"EPSG\",\"3395\"]]"; - } - } -}