From e8e7cacac071ea5055fa0947f610bce8afe02cd4 Mon Sep 17 00:00:00 2001 From: ClemensFischer Date: Sun, 25 Jan 2026 18:04:59 +0100 Subject: [PATCH] Abstract MapProjection.RelativeScale --- MapControl/Shared/MapProjection.cs | 10 +-- .../Shared/TransverseMercatorProjection.cs | 72 ++++++++++--------- .../TransverseMercatorProjectionSnyder.cs | 10 +-- MapProjections/Shared/ProjNetMapProjection.cs | 6 ++ 4 files changed, 54 insertions(+), 44 deletions(-) diff --git a/MapControl/Shared/MapProjection.cs b/MapControl/Shared/MapProjection.cs index d21a6abd..fc694c6d 100644 --- a/MapControl/Shared/MapProjection.cs +++ b/MapControl/Shared/MapProjection.cs @@ -108,7 +108,7 @@ namespace MapControl /// Gets the relative scale at the specified geographic coordinates. /// The returned Matrix represents the local distortion of the map projection. /// - public virtual Matrix RelativeScale(double latitude, double longitude) => new Matrix(1d, 0d, 0d, 1d, 0d, 0d); + public abstract Matrix RelativeScale(double latitude, double longitude); /// /// Transforms geographic coordinates to a Point in projected map coordinates. @@ -143,7 +143,7 @@ namespace MapControl /// Transforms a BoundingBox in geographic coordinates to a Rect in projected map coordinates. /// Returns null when the BoundingBox can not be transformed. /// - public virtual Rect? BoundingBoxToMap(BoundingBox boundingBox) + public Rect? BoundingBoxToMap(BoundingBox boundingBox) { var southWest = LocationToMap(boundingBox.South, boundingBox.West); var northEast = LocationToMap(boundingBox.North, boundingBox.East); @@ -153,9 +153,9 @@ namespace MapControl /// /// Transforms a Rect in projected map coordinates to a BoundingBox in geographic coordinates. - /// Returns null when the MapRect can not be transformed. + /// Returns null when the Rect can not be transformed. /// - public virtual BoundingBox MapToBoundingBox(Rect rect) + public BoundingBox MapToBoundingBox(Rect rect) { var southWest = MapToLocation(rect.X, rect.Y); var northEast = MapToLocation(rect.X + rect.Width, rect.Y + rect.Height); @@ -167,7 +167,7 @@ namespace MapControl /// Transforms a LatLonBox in geographic coordinates to a rotated Rect in projected map coordinates. /// Returns null when the LatLonBox can not be transformed. /// - public virtual (Rect?, double) LatLonBoxToMap(LatLonBox latLonBox) + public (Rect?, double) LatLonBoxToMap(LatLonBox latLonBox) { Rect? rect = null; var rotation = 0d; diff --git a/MapControl/Shared/TransverseMercatorProjection.cs b/MapControl/Shared/TransverseMercatorProjection.cs index 7654ff45..5ab6a681 100644 --- a/MapControl/Shared/TransverseMercatorProjection.cs +++ b/MapControl/Shared/TransverseMercatorProjection.cs @@ -1,6 +1,7 @@ using System; #if WPF using System.Windows; +using System.Windows.Media; #elif AVALONIA using Avalonia; #endif @@ -59,6 +60,42 @@ namespace MapControl Flattening = Wgs84Flattening; } + public override Matrix RelativeScale(double latitude, double longitude) + { + // h = k/k0 for relative scale + var h = 1d; + +#if TM_RELATIVE_SCALE // relative scale is usually < 1.001 and hence neglectable + // φ + 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 cosLambda = Math.Cos(dLambda); + // ξ' + var xi_ = Math.Atan2(t, cosLambda); + // η' + var eta_ = Atanh(Math.Sin(dLambda) / Math.Sqrt(1d + t * t)); + // σ + var sigma = 1 + + 2d * a1 * Math.Cos(2d * xi_) * Math.Cosh(2d * eta_) + + 4d * a2 * Math.Cos(4d * xi_) * Math.Cosh(4d * eta_) + + 6d * a3 * Math.Cos(6d * xi_) * Math.Cosh(6d * eta_); + // τ + var tau = + 2d * a1 * Math.Sin(2d * xi_) * Math.Sinh(2d * eta_) + + 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 u = (1d - n) / (1d + n) * Math.Tan(phi); + h = f1 * Math.Sqrt((1d + u * u) * (sigma * sigma + tau * tau) / (t * t + cosLambda * cosLambda)); +#endif + return new Matrix(h, 0d, 0d, h, 0d, 0d); + } + public override Point? LocationToMap(double latitude, double longitude) { // φ @@ -126,40 +163,5 @@ namespace MapControl #else private static double Atanh(double x) => Math.Atanh(x); #endif - /* - * Relative scale is usually < 1.001 and hence neglectable. - * - public override Matrix RelativeScale(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 cosLambda = Math.Cos(dLambda); - // ξ' - var xi_ = Math.Atan2(t, cosLambda); - // η' - var eta_ = Atanh(Math.Sin(dLambda) / Math.Sqrt(1d + t * t)); - // σ - var sigma = 1 + - 2d * a1 * Math.Cos(2d * xi_) * Math.Cosh(2d * eta_) + - 4d * a2 * Math.Cos(4d * xi_) * Math.Cosh(4d * eta_) + - 6d * a3 * Math.Cos(6d * xi_) * Math.Cosh(6d * eta_); - // τ - var tau = - 2d * a1 * Math.Sin(2d * xi_) * Math.Sinh(2d * eta_) + - 4d * a2 * Math.Sin(4d * xi_) * Math.Sinh(4d * eta_) + - 6d * a3 * Math.Sin(6d * xi_) * Math.Sinh(6d * eta_); - - // h = k/k0 for relative scale - var n = Flattening / (2d - Flattening); - var u = (1d - n) / (1d + n) * Math.Tan(phi); - var h = f1 * Math.Sqrt((1d + u * u) * (sigma * sigma + tau * tau) / (t * t + cosLambda * cosLambda)); - - return new Matrix(h, 0d, 0d, h, 0d, 0d); - }*/ } } diff --git a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs index 2dc924b7..c015dac0 100644 --- a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs +++ b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs @@ -47,8 +47,10 @@ namespace MapControl public override Matrix RelativeScale(double latitude, double longitude) { - var k = 1d; // omit k0 for relative scale + // h = k/k0 for relative scale + var h = 1d; +#if TM_RELATIVE_SCALE // relative scale is usually < 1.001 and hence neglectable if (latitude > -90d && latitude < 90d) { var phi = latitude * Math.PI / 180d; @@ -64,12 +66,12 @@ namespace MapControl var A4 = A2 * A2; var A6 = A2 * A4; - k *= 1d + (1d + C) * A2 / 2d + + h = 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) } - - return new Matrix(k, 0d, 0d, k, 0d, 0d); +#endif + return new Matrix(h, 0d, 0d, h, 0d, 0d); } public override Point? LocationToMap(double latitude, double longitude) diff --git a/MapProjections/Shared/ProjNetMapProjection.cs b/MapProjections/Shared/ProjNetMapProjection.cs index f3d3884c..f18c520c 100644 --- a/MapProjections/Shared/ProjNetMapProjection.cs +++ b/MapProjections/Shared/ProjNetMapProjection.cs @@ -3,6 +3,7 @@ using ProjNet.CoordinateSystems.Transformations; using System; #if WPF using System.Windows; +using System.Windows.Media; #elif AVALONIA using Avalonia; #endif @@ -70,6 +71,11 @@ namespace MapControl.Projections public MathTransform MapToLocationTransform { get; private set; } + public override Matrix RelativeScale(double latitude, double longitude) + { + return new Matrix(1d, 0d, 0d, 1d, 0d, 0d); + } + public override Point? LocationToMap(double latitude, double longitude) { if (LocationToMapTransform == null)