From 4ad9f2ea2a327c91b52f3a7b7889823b4eac6fc8 Mon Sep 17 00:00:00 2001 From: ClemensFischer Date: Thu, 29 Jan 2026 21:36:11 +0100 Subject: [PATCH] MapProjection.GridConvergence, rotation instead of Matrix --- MapControl/Shared/MapPanel.cs | 31 ++--- MapControl/Shared/MapProjection.cs | 46 ++++--- .../Shared/PolarStereographicProjection.cs | 58 +++++---- .../Shared/TransverseMercatorProjection.cs | 94 +++++++++----- .../TransverseMercatorProjectionSnyder.cs | 117 ++++++++++++------ MapProjections/Shared/Ed50UtmProjection.cs | 6 + MapProjections/Shared/Etrs89UtmProjection.cs | 17 +-- MapProjections/Shared/Nad27UtmProjection.cs | 17 +-- MapProjections/Shared/Nad83UtmProjection.cs | 17 +-- MapProjections/Shared/ProjNetMapProjection.cs | 24 +++- .../Shared/ProjNetMapProjectionFactory.cs | 20 +-- .../Shared/WebMercatorProjection.cs | 12 +- MapProjections/Shared/Wgs84UpsProjections.cs | 12 +- MapProjections/Shared/Wgs84UtmProjection.cs | 26 +--- .../Shared/WorldMercatorProjection.cs | 14 +-- 15 files changed, 272 insertions(+), 239 deletions(-) diff --git a/MapControl/Shared/MapPanel.cs b/MapControl/Shared/MapPanel.cs index 905ac450..fcdf93ce 100644 --- a/MapControl/Shared/MapPanel.cs +++ b/MapControl/Shared/MapPanel.cs @@ -280,7 +280,7 @@ namespace MapControl if (mapRect.HasValue) { - ArrangeElement(element, mapRect.Value, null); + ArrangeElement(element, mapRect.Value, 0d); } else { @@ -288,9 +288,9 @@ namespace MapControl if (boundingBox != null) { - (var rect, var transform) = parentMap.MapProjection.BoundingBoxToMap(boundingBox); + (var rect, var rotation) = parentMap.MapProjection.BoundingBoxToMap(boundingBox); - ArrangeElement(element, rect, transform); + ArrangeElement(element, rect, rotation); } else { @@ -300,7 +300,7 @@ namespace MapControl } } - private void ArrangeElement(FrameworkElement element, Rect mapRect, Matrix? transform) + private void ArrangeElement(FrameworkElement element, Rect mapRect, double rotation) { var viewRect = GetViewRect(mapRect); @@ -308,28 +308,15 @@ namespace MapControl element.Height = viewRect.Height; element.Arrange(viewRect); - if (parentMap.ViewTransform.Rotation != 0d) - { - var t = transform ?? new Matrix(1d, 0d, 0d, 1d, 0d, 0d); - t.Rotate(parentMap.ViewTransform.Rotation); - transform = t; - } + rotation += parentMap.ViewTransform.Rotation; - if (element.RenderTransform is MatrixTransform matrixTransform && - !matrixTransform.Matrix.IsIdentity) // not default RenderTransform in WPF/UWP/WinUI + if (element.RenderTransform is RotateTransform rotateTransform) { - if (transform.HasValue) - { - matrixTransform.Matrix = transform.Value; - } - else - { - element.ClearValue(RenderTransformProperty); - } + rotateTransform.Angle = rotation; } - else if (transform.HasValue) + else if (rotation != 0d) { - element.SetRenderTransform(new MatrixTransform { Matrix = transform.Value }, true); + element.SetRenderTransform(new RotateTransform { Angle = rotation }, true); } } diff --git a/MapControl/Shared/MapProjection.cs b/MapControl/Shared/MapProjection.cs index 4ae4001a..c218a13c 100644 --- a/MapControl/Shared/MapProjection.cs +++ b/MapControl/Shared/MapProjection.cs @@ -57,10 +57,17 @@ namespace MapControl /// /// Gets the relative transform at the specified geographic coordinates. - /// The returned Matrix represents the local distortion of the map projection. + /// The returned Matrix represents the local relative scale and rotation. /// public abstract Matrix RelativeTransform(double latitude, double longitude); + /// + /// Gets the grid convergence angle in degrees at the specified projected map coordinates. + /// Used for rotating the Rect resulting from BoundingBoxToMap in non-normal-cylindrical + /// projections, i.e. Transverse Mercator and Polar Stereographic. + /// + public virtual double GridConvergence(double x, double y) => 0d; + /// /// Transforms geographic coordinates to a Point in projected map coordinates. /// @@ -88,33 +95,42 @@ namespace MapControl /// /// Transforms a BoundingBox in geographic coordinates to a Rect in projected map coordinates - /// with an optional transform Matrix. + /// with an optional rotation angle in degrees. /// - public (Rect, Matrix?) BoundingBoxToMap(BoundingBox boundingBox) + public (Rect, double) BoundingBoxToMap(BoundingBox boundingBox) { Rect rect; - Matrix? transform = null; + var rotation = 0d; + var southWest = LocationToMap(boundingBox.South, boundingBox.West); + var northEast = LocationToMap(boundingBox.North, boundingBox.East); if (IsNormalCylindrical) { - var southWest = LocationToMap(boundingBox.South, boundingBox.West); - var northEast = LocationToMap(boundingBox.North, boundingBox.East); - rect = new Rect(southWest.X, southWest.Y, northEast.X - southWest.X, northEast.Y - southWest.Y); } else { - var latitude = (boundingBox.South + boundingBox.North) / 2d; - var longitude = (boundingBox.West + boundingBox.East) / 2d; - var center = LocationToMap(latitude, longitude); - var width = MeterPerDegree * (boundingBox.East - boundingBox.West) * Math.Cos(latitude * Math.PI / 180d); - var height = MeterPerDegree * (boundingBox.North - boundingBox.South); + var southEast = LocationToMap(boundingBox.South, boundingBox.East); + var northWest = LocationToMap(boundingBox.North, boundingBox.West); + var south = new Point((southWest.X + southEast.X) / 2d, (southWest.Y + southEast.Y) / 2d); + var north = new Point((northWest.X + northEast.X) / 2d, (northWest.Y + northEast.Y) / 2d); + var centerX = (south.X + north.X) / 2d; + var centerY = (south.Y + north.Y) / 2d; + var dx = north.X - south.X; + var dy = north.Y - south.Y; + var dxS = southEast.X - southWest.X; + var dyS = southEast.Y - southWest.Y; + var dxN = northEast.X - northWest.X; + var dyN = northEast.Y - northWest.Y; + var width = (Math.Sqrt(dxS * dxS + dyS * dyS) + Math.Sqrt(dxN * dxN + dyN * dyN)) / 2d; + var height = Math.Sqrt(dx * dx + dy * dy); - rect = new Rect(center.X - width / 2d, center.Y - height / 2d, width, height); - transform = RelativeTransform(latitude, longitude); + rect = new Rect(centerX - width / 2d, centerY - height / 2d, width, height); + + rotation = -GridConvergence(centerX, centerY); // invert direction for RotateTransform } - return (rect, transform); + return (rect, rotation); } /// diff --git a/MapControl/Shared/PolarStereographicProjection.cs b/MapControl/Shared/PolarStereographicProjection.cs index ac9e4dcd..1c436ea7 100644 --- a/MapControl/Shared/PolarStereographicProjection.cs +++ b/MapControl/Shared/PolarStereographicProjection.cs @@ -21,32 +21,6 @@ namespace MapControl public double FalseNorthing { get; set; } = 2e6; public Hemisphere Hemisphere { get; set; } - public static Matrix RelativeScale(Hemisphere hemisphere, double flattening, double scaleFactor, double latitude, double longitude) - { - var sign = hemisphere == Hemisphere.North ? 1d : -1d; - 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); - return transform; - } - - public override Matrix RelativeTransform(double latitude, double longitude) - { - return RelativeScale(Hemisphere, Flattening, ScaleFactor, latitude, longitude); - } - public override Point LocationToMap(double latitude, double longitude) { var sign = Hemisphere == Hemisphere.North ? 1d : -1d; @@ -68,6 +42,27 @@ namespace MapControl return new Point(x + FalseEasting, y + FalseNorthing); } + public override Matrix RelativeTransform(double latitude, double longitude) + { + var sign = Hemisphere == Hemisphere.North ? 1d : -1d; + 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); + return transform; + } + public override Location MapToLocation(double x, double y) { var sign = Hemisphere == Hemisphere.North ? 1d : -1d; @@ -84,6 +79,17 @@ namespace MapControl return new Location(sign * phi * 180d / Math.PI, sign * lambda * 180d / Math.PI); } + + public override double GridConvergence(double x, double y) + { + var sign = Hemisphere == Hemisphere.North ? 1d : -1d; + x = sign * (x - FalseEasting); + y = sign * (y - FalseNorthing); + + var lambda = Math.Atan2(x, -y); // p.162 (20-16) + + return lambda * 180d / Math.PI; + } } /// diff --git a/MapControl/Shared/TransverseMercatorProjection.cs b/MapControl/Shared/TransverseMercatorProjection.cs index 387967f0..c63f60ed 100644 --- a/MapControl/Shared/TransverseMercatorProjection.cs +++ b/MapControl/Shared/TransverseMercatorProjection.cs @@ -60,6 +60,35 @@ namespace MapControl Flattening = Wgs84Flattening; } + public override Point LocationToMap(double latitude, double longitude) + { + // φ + var phi = latitude * Math.PI / 180d; + var sinPhi = Math.Sin(phi); + // t + var t = Math.Sinh(Atanh(sinPhi) - f2 * Atanh(f2 * sinPhi)); + // λ - λ0 + var dLambda = (longitude - CentralMeridian) * Math.PI / 180d; + // ξ' + var xi_ = Math.Atan2(t, Math.Cos(dLambda)); + // η' + var eta_ = Atanh(Math.Sin(dLambda) / Math.Sqrt(1d + t * t)); + // k0 * A + var k0A = ScaleFactor * EquatorialRadius * f1; + + var x = FalseEasting + k0A * (eta_ + + a1 * Math.Cos(2d * xi_) * Math.Sinh(2d * eta_) + + a2 * Math.Cos(4d * xi_) * Math.Sinh(4d * eta_) + + a3 * Math.Cos(6d * xi_) * Math.Sinh(6d * eta_)); + + var y = FalseNorthing + k0A * (xi_ + + a1 * Math.Sin(2d * xi_) * Math.Cosh(2d * eta_) + + a2 * Math.Sin(4d * xi_) * Math.Cosh(4d * eta_) + + a3 * Math.Sin(6d * xi_) * Math.Cosh(6d * eta_)); + + return new Point(x, y); + } + public override Matrix RelativeTransform(double latitude, double longitude) { // φ @@ -100,35 +129,6 @@ namespace MapControl return transform; } - public override Point LocationToMap(double latitude, double longitude) - { - // φ - var phi = latitude * Math.PI / 180d; - var sinPhi = Math.Sin(phi); - // t - var t = Math.Sinh(Atanh(sinPhi) - f2 * Atanh(f2 * sinPhi)); - // λ - λ0 - var dLambda = (longitude - CentralMeridian) * Math.PI / 180d; - // ξ' - var xi_ = Math.Atan2(t, Math.Cos(dLambda)); - // η' - var eta_ = Atanh(Math.Sin(dLambda) / Math.Sqrt(1d + t * t)); - // k0 * A - var k0A = ScaleFactor * EquatorialRadius * f1; - - var x = FalseEasting + k0A * (eta_ + - a1 * Math.Cos(2d * xi_) * Math.Sinh(2d * eta_) + - a2 * Math.Cos(4d * xi_) * Math.Sinh(4d * eta_) + - a3 * Math.Cos(6d * xi_) * Math.Sinh(6d * eta_)); - - var y = FalseNorthing + k0A * (xi_ + - a1 * Math.Sin(2d * xi_) * Math.Cosh(2d * eta_) + - a2 * Math.Sin(4d * xi_) * Math.Cosh(4d * eta_) + - a3 * Math.Sin(6d * xi_) * Math.Cosh(6d * eta_)); - - return new Point(x, y); - } - public override Location MapToLocation(double x, double y) { // k0 * A @@ -162,6 +162,42 @@ namespace MapControl dLambda * 180d / Math.PI + CentralMeridian); } + public override double GridConvergence(double x, double y) + { + // k0 * A + var k0A = ScaleFactor * EquatorialRadius * f1; + // ξ + var xi = (y - FalseNorthing) / k0A; + // η + var eta = (x - FalseEasting) / k0A; + // ξ' + var xi_ = xi - + b1 * Math.Sin(2d * xi) * Math.Cosh(2d * eta) - + b2 * Math.Sin(4d * xi) * Math.Cosh(4d * eta) - + b3 * Math.Sin(6d * xi) * Math.Cosh(6d * eta); + // η' + var eta_ = eta - + b1 * Math.Cos(2d * xi) * Math.Sinh(2d * eta) - + b2 * Math.Cos(4d * xi) * Math.Sinh(4d * eta) - + b3 * Math.Cos(6d * xi) * Math.Sinh(6d * eta); + // σ' + var sigma_ = 1 - + 2d * b1 * Math.Cos(2d * xi) * Math.Cosh(2d * eta) + + 4d * b2 * Math.Cos(4d * xi) * Math.Cosh(4d * eta) + + 6d * b3 * Math.Cos(6d * xi) * Math.Cosh(6d * eta); + // τ' + var tau_ = + 2d * b1 * Math.Sin(2d * xi) * Math.Sinh(2d * eta) + + 4d * b2 * Math.Sin(4d * xi) * Math.Sinh(4d * eta) + + 6d * b3 * Math.Sin(6d * xi) * Math.Sinh(6d * eta); + + var tanXi_tanhEta_ = Math.Tan(xi_) * Math.Tanh(eta_); + // γ + var gamma = Math.Atan2(tau_ + sigma_ * tanXi_tanhEta_, sigma_ + -tau_ * tanXi_tanhEta_); + + return gamma * 180d / Math.PI; + } + #if NETFRAMEWORK private static double Atanh(double x) => Math.Log((1d + x) / (1d - x)) / 2d; #else diff --git a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs index 3bb06244..5cef2733 100644 --- a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs +++ b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs @@ -45,41 +45,6 @@ namespace MapControl e6 * 35d / 3072d * Math.Sin(6d * phi)); // (3-21) } - public override Matrix RelativeTransform(double latitude, double longitude) - { - var k = ScaleFactor; - var gamma = 0d; // γ, meridian convergence angle - - if (latitude > -90d && latitude < 90d) - { - var phi = latitude * Math.PI / 180d; - var sinPhi = Math.Sin(phi); - var cosPhi = Math.Cos(phi); - var tanPhi = sinPhi / cosPhi; - var dLambda = (longitude - CentralMeridian) * Math.PI / 180d; - - var e2 = (2d - Flattening) * Flattening; - var e_2 = e2 / (1d - e2); // (8-12) - var T = tanPhi * tanPhi; // (8-13) - var C = e_2 * cosPhi * cosPhi; // (8-14) - var A = dLambda * cosPhi; // (8-15) - var A2 = A * A; - var A4 = A2 * A2; - var A6 = A2 * A4; - - k *= 1d + (1d + C) * A2 / 2d + - (5d - 4d * T + 42d * C + 13d * C * C - 28d * e_2) * A4 / 24d + - (61d - 148d * T + 16 * T * T) * A6 / 720d; // (8-11) - - gamma = Math.Atan(Math.Tan(dLambda) * sinPhi); - } - - var transform = new Matrix(k, 0d, 0d, k, 0d, 0d); - transform.Rotate(-gamma * 180d / Math.PI); - - return transform; - } - public override Point LocationToMap(double latitude, double longitude) { var phi = latitude * Math.PI / 180d; @@ -118,6 +83,41 @@ namespace MapControl return new Point(x + FalseEasting, y + FalseNorthing); } + public override Matrix RelativeTransform(double latitude, double longitude) + { + var k = ScaleFactor; + var gamma = 0d; // γ, https://en.wikipedia.org/wiki/Transverse_Mercator_projection#Convergence + + if (latitude > -90d && latitude < 90d) + { + var phi = latitude * Math.PI / 180d; + var sinPhi = Math.Sin(phi); + var cosPhi = Math.Cos(phi); + var tanPhi = sinPhi / cosPhi; + var dLambda = (longitude - CentralMeridian) * Math.PI / 180d; + + var e2 = (2d - Flattening) * Flattening; + var e_2 = e2 / (1d - e2); // (8-12) + var T = tanPhi * tanPhi; // (8-13) + var C = e_2 * cosPhi * cosPhi; // (8-14) + var A = dLambda * cosPhi; // (8-15) + var A2 = A * A; + var A4 = A2 * A2; + var A6 = A2 * A4; + + k *= 1d + (1d + C) * A2 / 2d + + (5d - 4d * T + 42d * C + 13d * C * C - 28d * e_2) * A4 / 24d + + (61d - 148d * T + 16 * T * T) * A6 / 720d; // (8-11) + + gamma = Math.Atan(Math.Tan(dLambda) * sinPhi) * 180d / Math.PI; + } + + var transform = new Matrix(k, 0d, 0d, k, 0d, 0d); + transform.Rotate(-gamma); + + return transform; + } + public override Location MapToLocation(double x, double y) { var e2 = (2d - Flattening) * Flattening; @@ -164,5 +164,52 @@ namespace MapControl phi * 180d / Math.PI, dLambda * 180d / Math.PI + CentralMeridian); } + + public override double GridConvergence(double x, double y) + { + var e2 = (2d - Flattening) * Flattening; + var e4 = e2 * e2; + var e6 = e2 * e4; + var s = Math.Sqrt(1d - e2); + var e1 = (1d - s) / (1d + s); // (3-24) + var e12 = e1 * e1; + var e13 = e1 * e12; + var e14 = e1 * e13; + + var M = M0 + (y - FalseNorthing) / ScaleFactor; // (8-20) + var mu = M / (EquatorialRadius * (1d - e2 / 4d - e4 * 3d / 64d - e6 * 5d / 256d)); // (7-19) + var phi1 = mu + + (e1 * 3d / 2d - e13 * 27d / 32d) * Math.Sin(2d * mu) + + (e12 * 21d / 16d - e14 * 55d / 32d) * Math.Sin(4d * mu) + + e13 * 151d / 96d * Math.Sin(6d * mu) + + e14 * 1097d / 512d * Math.Sin(8d * mu); // (3-26) + + var sinPhi1 = Math.Sin(phi1); + var cosPhi1 = Math.Cos(phi1); + var tanPhi1 = sinPhi1 / cosPhi1; + + var e_2 = e2 / (1d - e2); // (8-12) + var C1 = e_2 * cosPhi1 * cosPhi1; // (8-21) + var T1 = sinPhi1 * sinPhi1 / (cosPhi1 * cosPhi1); // (8-22) + s = Math.Sqrt(1d - e2 * sinPhi1 * sinPhi1); + var N1 = EquatorialRadius / s; // (8-23) + var R1 = EquatorialRadius * (1d - e2) / (s * s * s); // (8-24) + var D = (x - FalseEasting) / (N1 * ScaleFactor); // (8-25) + var D2 = D * D; + var D3 = D * D2; + var D4 = D * D3; + var D5 = D * D4; + var D6 = D * D5; + + var phi = phi1 - N1 * tanPhi1 / R1 * (D2 / 2d - (5d + 3d * T1 + 10d * C1 - 4d * C1 * C1 - 9d * e_2) * D4 / 24d + + (61d + 90d * T1 + 45d * T1 * T1 + 298 * C1 - 3d * C1 * C1 - 252d * e_2) * D6 / 720d); // (8-17) + + var dLambda = (D - (1d + 2d * T1 + C1) * D3 / 6d + + (5d - 2d * C1 - 3d * C1 * C1 + 28d * T1 + 24d * T1 * T1 + 8d * e_2) * D5 / 120d) / cosPhi1; // (8-18) + + // γ, https://en.wikipedia.org/wiki/Transverse_Mercator_projection#Convergence + // + return Math.Atan(Math.Tan(dLambda) * Math.Sin(phi)) * 180d / Math.PI; + } } } diff --git a/MapProjections/Shared/Ed50UtmProjection.cs b/MapProjections/Shared/Ed50UtmProjection.cs index a11cb20a..0ef54a3f 100644 --- a/MapProjections/Shared/Ed50UtmProjection.cs +++ b/MapProjections/Shared/Ed50UtmProjection.cs @@ -15,6 +15,12 @@ namespace MapControl.Projections public int Zone { get; } public Ed50UtmProjection(int zone) + : base(new TransverseMercatorProjection + { + EquatorialRadius = 6378388d, + Flattening = 297, + CentralMeridian = 6 * zone - 183, + }) { if (zone < FirstZone || zone > LastZone) { diff --git a/MapProjections/Shared/Etrs89UtmProjection.cs b/MapProjections/Shared/Etrs89UtmProjection.cs index 7e281f21..e1b21514 100644 --- a/MapProjections/Shared/Etrs89UtmProjection.cs +++ b/MapProjections/Shared/Etrs89UtmProjection.cs @@ -1,26 +1,15 @@ -using System; - -namespace MapControl.Projections +namespace MapControl.Projections { /// - /// ETRS89 Universal Transverse Mercator Projection. + /// ETRS89 Universal Transverse Mercator Projection - EPSG:25828 to EPSG:25838. /// public class Etrs89UtmProjection : ProjNetMapProjection { - public const int FirstZone = 28; - public const int LastZone = 38; - public const int FirstZoneEpsgCode = 25800 + FirstZone; - public const int LastZoneEpsgCode = 25800 + LastZone; - public int Zone { get; } public Etrs89UtmProjection(int zone) + : base(new MapControl.Etrs89UtmProjection(zone)) { - if (zone < FirstZone || zone > LastZone) - { - throw new ArgumentException($"Invalid ETRS89 UTM zone {zone}.", nameof(zone)); - } - Zone = zone; CoordinateSystemWkt = $"PROJCS[\"ETRS89 / UTM zone {zone}N\"," + diff --git a/MapProjections/Shared/Nad27UtmProjection.cs b/MapProjections/Shared/Nad27UtmProjection.cs index 080583a6..098b5e42 100644 --- a/MapProjections/Shared/Nad27UtmProjection.cs +++ b/MapProjections/Shared/Nad27UtmProjection.cs @@ -1,26 +1,15 @@ -using System; - -namespace MapControl.Projections +namespace MapControl.Projections { /// - /// NAD27 Universal Transverse Mercator Projection. + /// NAD27 Universal Transverse Mercator Projection - EPSG:26701 to EPSG:26722. /// public class Nad27UtmProjection : ProjNetMapProjection { - public const int FirstZone = 1; - public const int LastZone = 22; - public const int FirstZoneEpsgCode = 26700 + FirstZone; - public const int LastZoneEpsgCode = 26700 + LastZone; - public int Zone { get; } public Nad27UtmProjection(int zone) + : base(new MapControl.Nad27UtmProjection(zone)) { - if (zone < FirstZone || zone > LastZone) - { - throw new ArgumentException($"Invalid NAD27 UTM zone {zone}.", nameof(zone)); - } - Zone = zone; CoordinateSystemWkt = $"PROJCS[\"NAD27 / UTM zone {zone}N\"," + diff --git a/MapProjections/Shared/Nad83UtmProjection.cs b/MapProjections/Shared/Nad83UtmProjection.cs index e2cb20f1..9c39dccf 100644 --- a/MapProjections/Shared/Nad83UtmProjection.cs +++ b/MapProjections/Shared/Nad83UtmProjection.cs @@ -1,26 +1,15 @@ -using System; - -namespace MapControl.Projections +namespace MapControl.Projections { /// - /// NAD83 Universal Transverse Mercator Projection. + /// NAD83 Universal Transverse Mercator Projection - EPSG:26901 to EPSG:26923. /// public class Nad83UtmProjection : ProjNetMapProjection { - public const int FirstZone = 1; - public const int LastZone = 23; - public const int FirstZoneEpsgCode = 26900 + FirstZone; - public const int LastZoneEpsgCode = 26900 + LastZone; - public int Zone { get; } public Nad83UtmProjection(int zone) + : base(new MapControl.Nad83UtmProjection(zone)) { - if (zone < FirstZone || zone > LastZone) - { - throw new ArgumentException($"Invalid NAD83 UTM zone {zone}.", nameof(zone)); - } - Zone = zone; CoordinateSystemWkt = $"PROJCS[\"NAD83 / UTM zone {zone}N\"," + diff --git a/MapProjections/Shared/ProjNetMapProjection.cs b/MapProjections/Shared/ProjNetMapProjection.cs index 762b2c97..9c813939 100644 --- a/MapProjections/Shared/ProjNetMapProjection.cs +++ b/MapProjections/Shared/ProjNetMapProjection.cs @@ -15,8 +15,11 @@ namespace MapControl.Projections /// public class ProjNetMapProjection : MapProjection { - protected ProjNetMapProjection() + protected MapProjection FallbackProjection { get; } + + protected ProjNetMapProjection(MapProjection fallbackProjection) { + FallbackProjection = fallbackProjection; } public ProjNetMapProjection(string coordinateSystemWkt) @@ -71,11 +74,6 @@ namespace MapControl.Projections public MathTransform MapToLocationTransform { get; private set; } - public override Matrix RelativeTransform(double latitude, double longitude) - { - return new Matrix(1d, 0d, 0d, 1d, 0d, 0d); - } - public override Point LocationToMap(double latitude, double longitude) { if (LocationToMapTransform == null) @@ -99,5 +97,19 @@ namespace MapControl.Projections return new Location(coordinate[1], coordinate[0]); } + + public override Matrix RelativeTransform(double latitude, double longitude) + { + return FallbackProjection != null + ? FallbackProjection.RelativeTransform(latitude, longitude) + : new Matrix(1d, 0d, 0d, 1d, 0d, 0d); + } + + public override double GridConvergence(double x, double y) + { + return FallbackProjection != null + ? FallbackProjection.GridConvergence(x, y) + : 0d; + } } } diff --git a/MapProjections/Shared/ProjNetMapProjectionFactory.cs b/MapProjections/Shared/ProjNetMapProjectionFactory.cs index 7537edc5..6964ddee 100644 --- a/MapProjections/Shared/ProjNetMapProjectionFactory.cs +++ b/MapProjections/Shared/ProjNetMapProjectionFactory.cs @@ -47,16 +47,16 @@ namespace MapControl.Projections { var c when c is >= Ed50UtmProjection.FirstZoneEpsgCode and <= Ed50UtmProjection.LastZoneEpsgCode => new Ed50UtmProjection(c % 100), - var c when c is >= Etrs89UtmProjection.FirstZoneEpsgCode - and <= Etrs89UtmProjection.LastZoneEpsgCode => new Etrs89UtmProjection(c % 100), - var c when c is >= Nad27UtmProjection.FirstZoneEpsgCode - and <= Nad27UtmProjection.LastZoneEpsgCode => new Nad27UtmProjection(c % 100), - 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), - var c when c is >= Wgs84UtmProjection.FirstZoneSouthEpsgCode - and <= Wgs84UtmProjection.LastZoneSouthEpsgCode => new Wgs84UtmProjection(c % 100, Hemisphere.South), + var c when c is >= MapControl.Etrs89UtmProjection.FirstZoneEpsgCode + and <= MapControl.Etrs89UtmProjection.LastZoneEpsgCode => new Etrs89UtmProjection(c % 100), + var c when c is >= MapControl.Nad27UtmProjection.FirstZoneEpsgCode + and <= MapControl.Nad27UtmProjection.LastZoneEpsgCode => new Nad27UtmProjection(c % 100), + 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), + var c when c is >= MapControl.Wgs84UtmProjection.FirstZoneSouthEpsgCode + and <= MapControl.Wgs84UtmProjection.LastZoneSouthEpsgCode => new Wgs84UtmProjection(c % 100, Hemisphere.South), _ => base.GetProjection(epsgCode) }; } diff --git a/MapProjections/Shared/WebMercatorProjection.cs b/MapProjections/Shared/WebMercatorProjection.cs index 9c32b71d..37659a80 100644 --- a/MapProjections/Shared/WebMercatorProjection.cs +++ b/MapProjections/Shared/WebMercatorProjection.cs @@ -1,8 +1,4 @@ using ProjNet.CoordinateSystems; -using System; -#if WPF -using System.Windows.Media; -#endif namespace MapControl.Projections { @@ -13,15 +9,9 @@ namespace MapControl.Projections public class WebMercatorProjection : ProjNetMapProjection { public WebMercatorProjection() + : base(new MapControl.WebMercatorProjection()) { CoordinateSystem = ProjectedCoordinateSystem.WebMercator; } - - public override Matrix RelativeTransform(double latitude, double longitude) - { - var k = 1d / Math.Cos(latitude * Math.PI / 180d); // p.44 (7-3) - - return new Matrix(k, 0d, 0d, k, 0d, 0d); - } } } diff --git a/MapProjections/Shared/Wgs84UpsProjections.cs b/MapProjections/Shared/Wgs84UpsProjections.cs index 055972d9..e72da80e 100644 --- a/MapProjections/Shared/Wgs84UpsProjections.cs +++ b/MapProjections/Shared/Wgs84UpsProjections.cs @@ -7,6 +7,7 @@ namespace MapControl.Projections public class Wgs84UpsNorthProjection : ProjNetMapProjection { public Wgs84UpsNorthProjection() + : base(new MapControl.Wgs84UpsNorthProjection()) { CoordinateSystemWkt = "PROJCS[\"WGS 84 / UPS North (N,E)\"," + @@ -20,16 +21,12 @@ namespace MapControl.Projections "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + "AUTHORITY[\"EPSG\",\"32661\"]]"; } - - public override Matrix RelativeTransform(double latitude, double longitude) - { - return PolarStereographicProjection.RelativeScale(Hemisphere.North, Wgs84Flattening, 0.994, latitude, longitude); - } } public class Wgs84UpsSouthProjection : ProjNetMapProjection { public Wgs84UpsSouthProjection() + : base(new MapControl.Wgs84UpsSouthProjection()) { CoordinateSystemWkt = "PROJCS[\"WGS 84 / UPS South (N,E)\"," + @@ -43,10 +40,5 @@ namespace MapControl.Projections "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," + "AUTHORITY[\"EPSG\",\"32761\"]]"; } - - public override Matrix RelativeTransform(double latitude, double longitude) - { - return PolarStereographicProjection.RelativeScale(Hemisphere.North, Wgs84Flattening, 0.994, latitude, longitude); - } } } diff --git a/MapProjections/Shared/Wgs84UtmProjection.cs b/MapProjections/Shared/Wgs84UtmProjection.cs index d0d2db57..b6da1918 100644 --- a/MapProjections/Shared/Wgs84UtmProjection.cs +++ b/MapProjections/Shared/Wgs84UtmProjection.cs @@ -1,35 +1,19 @@ using ProjNet.CoordinateSystems; -using System; namespace MapControl.Projections { /// - /// WGS84 Universal Transverse Mercator Projection. + /// WGS84 Universal Transverse Mercator Projection - + /// EPSG:32601 to EPSG:32660 and EPSG:32701 to EPSG:32760. /// public class Wgs84UtmProjection : ProjNetMapProjection { - public const int FirstZone = 1; - public const int LastZone = 60; - public const int FirstZoneNorthEpsgCode = 32600 + FirstZone; - public const int LastZoneNorthEpsgCode = 32600 + LastZone; - public const int FirstZoneSouthEpsgCode = 32700 + FirstZone; - public const int LastZoneSouthEpsgCode = 32700 + LastZone; - - public int Zone { get; private set; } - public Hemisphere Hemisphere { get; private set; } + public int Zone { get; } + public Hemisphere Hemisphere { get; } public Wgs84UtmProjection(int zone, Hemisphere hemisphere) + : base(new MapControl.Wgs84UtmProjection(zone, hemisphere)) { - SetZone(zone, hemisphere); - } - - protected void SetZone(int zone, Hemisphere hemisphere) - { - if (zone < FirstZone || zone > LastZone) - { - throw new ArgumentException($"Invalid WGS84 UTM zone {zone}.", nameof(zone)); - } - Zone = zone; Hemisphere = hemisphere; CoordinateSystem = ProjectedCoordinateSystem.WGS84_UTM(zone, hemisphere == Hemisphere.North); diff --git a/MapProjections/Shared/WorldMercatorProjection.cs b/MapProjections/Shared/WorldMercatorProjection.cs index 63c9708c..77a512ca 100644 --- a/MapProjections/Shared/WorldMercatorProjection.cs +++ b/MapProjections/Shared/WorldMercatorProjection.cs @@ -1,8 +1,4 @@ -#if WPF -using System.Windows.Media; -#endif - -namespace MapControl.Projections +namespace MapControl.Projections { /// /// Elliptical Mercator Projection implemented by setting the WKT property of a ProjNetMapProjection. @@ -11,6 +7,7 @@ namespace MapControl.Projections public class WorldMercatorProjection : ProjNetMapProjection { public WorldMercatorProjection() + : base(new MapControl.WorldMercatorProjection()) { CoordinateSystemWkt = "PROJCS[\"WGS 84 / World Mercator\"," + @@ -26,12 +23,5 @@ namespace MapControl.Projections "AXIS[\"Northing\",NORTH]," + "AUTHORITY[\"EPSG\",\"3395\"]]"; } - - public override Matrix RelativeTransform(double latitude, double longitude) - { - var k = MapControl.WorldMercatorProjection.RelativeScale(latitude); - - return new Matrix(k, 0d, 0d, k, 0d, 0d); - } } }