From 3922297e7b78489a02cc0d649b3ba34c2e86774a Mon Sep 17 00:00:00 2001 From: ClemensFischer Date: Wed, 28 Jan 2026 09:52:11 +0100 Subject: [PATCH] Fixed RelativeTransform of azimuthal projections --- .../Shared/AzimuthalEquidistantProjection.cs | 21 ++++++++-------- MapControl/Shared/AzimuthalProjection.cs | 23 +++++++++++------ MapControl/Shared/GnomonicProjection.cs | 3 ++- MapControl/Shared/MapProjection.cs | 25 +++++++++++++++++++ MapControl/Shared/OrthographicProjection.cs | 3 ++- .../Shared/PolarStereographicProjection.cs | 8 +++--- MapControl/Shared/StereographicProjection.cs | 2 +- .../Shared/TransverseMercatorProjection.cs | 7 +++--- .../TransverseMercatorProjectionSnyder.cs | 10 +++----- .../Shared/Wgs84OrthographicProjection.cs | 3 ++- .../Shared/Wgs84StereographicProjection.cs | 2 +- MapProjections/Shared/Wgs84UpsProjections.cs | 4 +-- 12 files changed, 72 insertions(+), 39 deletions(-) diff --git a/MapControl/Shared/AzimuthalEquidistantProjection.cs b/MapControl/Shared/AzimuthalEquidistantProjection.cs index 193795b4..977382ca 100644 --- a/MapControl/Shared/AzimuthalEquidistantProjection.cs +++ b/MapControl/Shared/AzimuthalEquidistantProjection.cs @@ -29,21 +29,22 @@ namespace MapControl public override Matrix RelativeTransform(double latitude, double longitude) { var p = GetProjectedPoint(latitude, longitude); - - if (p.CosC == 1d) - { - return new Matrix(1d, 0d, 0d, 1d, 0d, 0d); - } + var scaleX = 1d; + var scaleY = 1d; if (p.CosC == -1d) { - return new Matrix(1d, 0d, 0d, double.PositiveInfinity, 0d, 0d); + scaleY = double.PositiveInfinity; + } + else if (p.CosC != 1d) + { + var c = Math.Acos(p.CosC); + var k = c / Math.Sin(c); // p.195 (25-2) + + (scaleX, scaleY) = p.RelativeScale(1d, k); } - var c = Math.Acos(p.CosC); - var k = c / Math.Sin(c); // p.195 (25-2) - - return p.RelativeScale(1d, k); + return RelativeTransform(latitude, longitude, scaleX, scaleY); } public override Point? LocationToMap(double latitude, double longitude) diff --git a/MapControl/Shared/AzimuthalProjection.cs b/MapControl/Shared/AzimuthalProjection.cs index 6acdb99c..9e54aa42 100644 --- a/MapControl/Shared/AzimuthalProjection.cs +++ b/MapControl/Shared/AzimuthalProjection.cs @@ -1,7 +1,4 @@ using System; -#if WPF -using System.Windows.Media; -#endif namespace MapControl { @@ -36,16 +33,26 @@ namespace MapControl CosC = Math.Min(Math.Max(cosC, -1d), 1d); // protect against rounding errors } - public Matrix RelativeScale(double radialScale, double perpendicularScale) + public (double, double) RelativeScale(double radialScale, double perpendicularScale) { - var scale = new Matrix(radialScale, 0d, 0d, perpendicularScale, 0d, 0d); + var scaleX = radialScale; + var scaleY = perpendicularScale; - if (radialScale != perpendicularScale) + if (scaleX != scaleY) { - scale.Rotate(-Math.Atan2(Y, X) * 180d / Math.PI); + var s = Math.Sqrt(X * X + Y * Y); + // sin and cos of azimuth from projection center, i.e. Atan2(Y/X) + var cos = X / s; + var sin = Y / s; + var x1 = scaleX * cos; + var y1 = scaleY * sin; + var x2 = scaleX * sin; + var y2 = scaleY * cos; + scaleX = Math.Sqrt(x1 * x1 + y1 * y1); + scaleY = Math.Sqrt(x2 * x2 + y2 * y2); } - return scale; + return (scaleX, scaleY); } } diff --git a/MapControl/Shared/GnomonicProjection.cs b/MapControl/Shared/GnomonicProjection.cs index 59d57c7c..c7823425 100644 --- a/MapControl/Shared/GnomonicProjection.cs +++ b/MapControl/Shared/GnomonicProjection.cs @@ -31,8 +31,9 @@ namespace MapControl var p = GetProjectedPoint(latitude, longitude); var k = 1d / p.CosC; // p.165 (22-3) var h = k * k; // p.165 (22-2) + (var scaleX, var scaleY) = p.RelativeScale(h, k); - return p.RelativeScale(h, k); + return RelativeTransform(latitude, longitude, scaleX, scaleY); } public override Point? LocationToMap(double latitude, double longitude) diff --git a/MapControl/Shared/MapProjection.cs b/MapControl/Shared/MapProjection.cs index 3aaa3fa8..57916e9c 100644 --- a/MapControl/Shared/MapProjection.cs +++ b/MapControl/Shared/MapProjection.cs @@ -207,5 +207,30 @@ namespace MapControl { return CrsId; } + + /// + /// Used by azimuthal projections, where local skewing can not be directly calculated. + /// + protected Matrix RelativeTransform(double latitude, double longitude, double scaleX, double scaleY) + { + var north = LocationToMap(latitude - 1e-3, longitude); + var south = LocationToMap(latitude + 1e-3, longitude); + var west = LocationToMap(latitude, longitude - 1e-3); + var east = LocationToMap(latitude, longitude + 1e-3); + var tanSkewX = 0d; + var tanSkewY = 0d; + + if (north.HasValue && south.HasValue && west.HasValue && east.HasValue) + { + var dx1 = east.Value.X - west.Value.X; + var dy1 = east.Value.Y - west.Value.Y; + var dx2 = north.Value.X - south.Value.X; + var dy2 = north.Value.Y - south.Value.Y; + tanSkewX = -dx2 / dy2; + tanSkewY = -dy1 / dx1; + } + + return new Matrix(scaleX, scaleX * tanSkewY, scaleY * tanSkewX, scaleY, 0d, 0d); + } } } diff --git a/MapControl/Shared/OrthographicProjection.cs b/MapControl/Shared/OrthographicProjection.cs index 263576a0..7750f285 100644 --- a/MapControl/Shared/OrthographicProjection.cs +++ b/MapControl/Shared/OrthographicProjection.cs @@ -29,8 +29,9 @@ namespace MapControl public override Matrix RelativeTransform(double latitude, double longitude) { var p = GetProjectedPoint(latitude, longitude); + (var scaleX, var scaleY) = p.RelativeScale(p.CosC, 1d); // p.149 (20-5), k == 1 - return p.RelativeScale(p.CosC, 1d); // p.149 (20-5), k == 1 + return RelativeTransform(latitude, longitude, scaleX, scaleY); } public override Point? LocationToMap(double latitude, double longitude) diff --git a/MapControl/Shared/PolarStereographicProjection.cs b/MapControl/Shared/PolarStereographicProjection.cs index 787b5d73..0e7fb09e 100644 --- a/MapControl/Shared/PolarStereographicProjection.cs +++ b/MapControl/Shared/PolarStereographicProjection.cs @@ -21,7 +21,7 @@ namespace MapControl public double FalseNorthing { get; set; } = 2e6; public Hemisphere Hemisphere { get; set; } - public static double RelativeScale(Hemisphere hemisphere, double flattening, double latitude) + public static double RelativeScale(Hemisphere hemisphere, double flattening, double scaleFactor, double latitude) { var sign = hemisphere == Hemisphere.North ? 1d : -1d; var phi = sign * latitude * Math.PI / 180d; @@ -32,8 +32,8 @@ namespace MapControl var t = Math.Tan(Math.PI / 4d - phi / 2d) / Math.Pow((1d - eSinPhi) / (1d + eSinPhi), e / 2d); // p.161 (15-9) - // r == ρ/(a*k0), omit k0 for relative scale - var r = 2d * t / Math.Sqrt(Math.Pow(1d + e, 1d + e) * Math.Pow(1d - e, 1d - e)); // p.161 (21-33) + // r == ρ/a + var r = scaleFactor * 2d * 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) return r / m; // p.161 (21-32) @@ -41,7 +41,7 @@ namespace MapControl public override Matrix RelativeTransform(double latitude, double longitude) { - var k = RelativeScale(Hemisphere, Flattening, latitude); + var k = RelativeScale(Hemisphere, Flattening, ScaleFactor, latitude); return new Matrix(k, 0d, 0d, k, 0d, 0d); } diff --git a/MapControl/Shared/StereographicProjection.cs b/MapControl/Shared/StereographicProjection.cs index 9b891793..84d0065f 100644 --- a/MapControl/Shared/StereographicProjection.cs +++ b/MapControl/Shared/StereographicProjection.cs @@ -31,7 +31,7 @@ namespace MapControl var p = GetProjectedPoint(latitude, longitude); var k = 2d / (1d + p.CosC); // p.157 (21-4), k0 == 1 - return p.RelativeScale(k, k); + return RelativeTransform(latitude, longitude, k, k); } public override Point? LocationToMap(double latitude, double longitude) diff --git a/MapControl/Shared/TransverseMercatorProjection.cs b/MapControl/Shared/TransverseMercatorProjection.cs index 3ab4da9c..70f12fc3 100644 --- a/MapControl/Shared/TransverseMercatorProjection.cs +++ b/MapControl/Shared/TransverseMercatorProjection.cs @@ -68,6 +68,7 @@ namespace MapControl // λ - λ0 var dLambda = (longitude - CentralMeridian) * Math.PI / 180d; var cosLambda = Math.Cos(dLambda); + var tanLambda = Math.Tan(dLambda); // t var t = Math.Sinh(Atanh(sinPhi) - f2 * Atanh(f2 * sinPhi)); var u = Math.Sqrt(1d + t * t); @@ -88,14 +89,12 @@ namespace MapControl var n = Flattening / (2d - Flattening); var m = (1d - n) / (1d + n) * Math.Tan(phi); - // h = k/k0 for relative scale - var h = f1 * Math.Sqrt((1d + m * m) * (sigma * sigma + tau * tau) / (t * t + cosLambda * cosLambda)); + var k = ScaleFactor * f1 * Math.Sqrt((1d + m * m) * (sigma * sigma + tau * tau) / (t * t + cosLambda * cosLambda)); - var tanLambda = Math.Tan(dLambda); // γ, meridian convergence angle var gamma = Math.Atan2(tau * u + sigma * t * tanLambda, sigma * u - tau * t * tanLambda); - var transform = new Matrix(h, 0d, 0d, h, 0d, 0d); + var transform = new Matrix(k, 0d, 0d, k, 0d, 0d); transform.Rotate(-gamma * 180d / Math.PI); return transform; diff --git a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs index e9ab561d..07bfd337 100644 --- a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs +++ b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs @@ -47,10 +47,8 @@ namespace MapControl public override Matrix RelativeTransform(double latitude, double longitude) { - // h = k/k0 for relative scale - var h = 1d; - // γ, meridian convergence angle - var gamma = 0d; + var k = ScaleFactor; + var gamma = 0d; // γ, meridian convergence angle if (latitude > -90d && latitude < 90d) { @@ -69,14 +67,14 @@ namespace MapControl var A4 = A2 * A2; var A6 = A2 * A4; - h = 1d + (1d + C) * A2 / 2d + + 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(h, 0d, 0d, h, 0d, 0d); + var transform = new Matrix(k, 0d, 0d, k, 0d, 0d); transform.Rotate(-gamma * 180d / Math.PI); return transform; diff --git a/MapProjections/Shared/Wgs84OrthographicProjection.cs b/MapProjections/Shared/Wgs84OrthographicProjection.cs index 400f6f60..f4a29487 100644 --- a/MapProjections/Shared/Wgs84OrthographicProjection.cs +++ b/MapProjections/Shared/Wgs84OrthographicProjection.cs @@ -37,8 +37,9 @@ namespace MapControl.Projections public override Matrix RelativeTransform(double latitude, double longitude) { var p = new AzimuthalProjection.ProjectedPoint(Center.Latitude, Center.Longitude, latitude, longitude); + (var scaleX, var scaleY) = p.RelativeScale(p.CosC, 1d); // p.149 (20-5), k == 1 - return p.RelativeScale(p.CosC, 1d); // p.149 (20-5), k == 1 + return RelativeTransform(latitude, longitude, scaleX, scaleY); } } } diff --git a/MapProjections/Shared/Wgs84StereographicProjection.cs b/MapProjections/Shared/Wgs84StereographicProjection.cs index c697007f..55cc7db9 100644 --- a/MapProjections/Shared/Wgs84StereographicProjection.cs +++ b/MapProjections/Shared/Wgs84StereographicProjection.cs @@ -39,7 +39,7 @@ namespace MapControl.Projections var p = new AzimuthalProjection.ProjectedPoint(Center.Latitude, Center.Longitude, latitude, longitude); var k = 2d / (1d + p.CosC); // p.157 (21-4), k0 == 1 - return p.RelativeScale(k, k); + return RelativeTransform(latitude, longitude, k, k); } } } diff --git a/MapProjections/Shared/Wgs84UpsProjections.cs b/MapProjections/Shared/Wgs84UpsProjections.cs index fdb6acb5..9ee1092d 100644 --- a/MapProjections/Shared/Wgs84UpsProjections.cs +++ b/MapProjections/Shared/Wgs84UpsProjections.cs @@ -23,7 +23,7 @@ namespace MapControl.Projections public override Matrix RelativeTransform(double latitude, double longitude) { - var k = PolarStereographicProjection.RelativeScale(Hemisphere.North, Wgs84Flattening, latitude); + var k = PolarStereographicProjection.RelativeScale(Hemisphere.North, Wgs84Flattening, 0.994, latitude); return new Matrix(k, 0d, 0d, k, 0d, 0d); } @@ -48,7 +48,7 @@ namespace MapControl.Projections public override Matrix RelativeTransform(double latitude, double longitude) { - var k = PolarStereographicProjection.RelativeScale(Hemisphere.South, Wgs84Flattening, latitude); + var k = PolarStereographicProjection.RelativeScale(Hemisphere.South, Wgs84Flattening, 0.994, latitude); return new Matrix(k, 0d, 0d, k, 0d, 0d); }