From 23a8e49efb87e908bdf4a9fe3fd85c6da9be34a4 Mon Sep 17 00:00:00 2001 From: ClemensFischer Date: Fri, 9 Jan 2026 08:13:07 +0100 Subject: [PATCH] Use mean radius in distance calculations --- MapControl/Shared/AzimuthalEquidistantProjection.cs | 4 ++-- MapControl/Shared/GnomonicProjection.cs | 4 ++-- MapControl/Shared/Location.cs | 4 ++-- MapControl/Shared/MapProjection.cs | 6 +++++- MapControl/Shared/OrthographicProjection.cs | 8 ++++---- MapControl/Shared/StereographicProjection.cs | 4 ++-- MapControl/Shared/WorldMercatorProjection.cs | 2 ++ MapProjections/Shared/WorldMercatorProjection.cs | 2 +- 8 files changed, 20 insertions(+), 14 deletions(-) diff --git a/MapControl/Shared/AzimuthalEquidistantProjection.cs b/MapControl/Shared/AzimuthalEquidistantProjection.cs index 5a1eb738..52027047 100644 --- a/MapControl/Shared/AzimuthalEquidistantProjection.cs +++ b/MapControl/Shared/AzimuthalEquidistantProjection.cs @@ -36,7 +36,7 @@ namespace MapControl Center.GetAzimuthDistance(latitude, longitude, out double azimuth, out double distance); - var mapDistance = distance * Wgs84EquatorialRadius; + var mapDistance = distance * Wgs84MeanRadius; return new Point(mapDistance * Math.Sin(azimuth), mapDistance * Math.Cos(azimuth)); } @@ -50,7 +50,7 @@ namespace MapControl var azimuth = Math.Atan2(x, y); var mapDistance = Math.Sqrt(x * x + y * y); - var distance = mapDistance / Wgs84EquatorialRadius; + var distance = mapDistance / Wgs84MeanRadius; return Center.GetLocation(azimuth, distance); } diff --git a/MapControl/Shared/GnomonicProjection.cs b/MapControl/Shared/GnomonicProjection.cs index 44699cbc..5013e35c 100644 --- a/MapControl/Shared/GnomonicProjection.cs +++ b/MapControl/Shared/GnomonicProjection.cs @@ -37,7 +37,7 @@ namespace MapControl Center.GetAzimuthDistance(latitude, longitude, out double azimuth, out double distance); var mapDistance = distance < Math.PI / 2d - ? Math.Tan(distance) * Wgs84EquatorialRadius + ? Math.Tan(distance) * Wgs84MeanRadius : double.PositiveInfinity; return new Point(mapDistance * Math.Sin(azimuth), mapDistance * Math.Cos(azimuth)); @@ -52,7 +52,7 @@ namespace MapControl var azimuth = Math.Atan2(x, y); var mapDistance = Math.Sqrt(x * x + y * y); - var distance = Math.Atan(mapDistance / Wgs84EquatorialRadius); + var distance = Math.Atan(mapDistance / Wgs84MeanRadius); return Center.GetLocation(azimuth, distance); } diff --git a/MapControl/Shared/Location.cs b/MapControl/Shared/Location.cs index 143db6ab..74f1a47a 100644 --- a/MapControl/Shared/Location.cs +++ b/MapControl/Shared/Location.cs @@ -115,7 +115,7 @@ namespace MapControl /// /// Calculates the great circle distance in meters between this and the specified Location. /// - public double GetDistance(Location location, double earthRadius = MapProjection.Wgs84EquatorialRadius) + public double GetDistance(Location location, double earthRadius = MapProjection.Wgs84MeanRadius) { GetAzimuthDistance(location.Latitude, location.Longitude, out _, out double distance); @@ -144,7 +144,7 @@ namespace MapControl /// /// Calculates the Location on a great circle at the specified azimuth in degrees and distance in meters from this Location. /// - public Location GetLocation(double azimuth, double distance, double earthRadius = MapProjection.Wgs84EquatorialRadius) + public Location GetLocation(double azimuth, double distance, double earthRadius = MapProjection.Wgs84MeanRadius) { return GetLocation(azimuth * Math.PI / 180d, distance / earthRadius); } diff --git a/MapControl/Shared/MapProjection.cs b/MapControl/Shared/MapProjection.cs index 8cb106bb..1bd8856e 100644 --- a/MapControl/Shared/MapProjection.cs +++ b/MapControl/Shared/MapProjection.cs @@ -29,7 +29,11 @@ namespace MapControl public const double Wgs84EquatorialRadius = 6378137d; public const double Wgs84MeterPerDegree = Wgs84EquatorialRadius * Math.PI / 180d; public const double Wgs84Flattening = 1d / 298.257223563; - public static readonly double Wgs84Eccentricity = Math.Sqrt((2d - Wgs84Flattening) * Wgs84Flattening); + + // Arithmetic mean radius (2*a + b) / 3 == (1 - f/3) * a + // https://en.wikipedia.org/wiki/Earth_radius#Arithmetic_mean_radius + // + public const double Wgs84MeanRadius = (1d - Wgs84Flattening / 3d) * Wgs84EquatorialRadius; public static MapProjectionFactory Factory { diff --git a/MapControl/Shared/OrthographicProjection.cs b/MapControl/Shared/OrthographicProjection.cs index 4623fcf4..3a764ef5 100644 --- a/MapControl/Shared/OrthographicProjection.cs +++ b/MapControl/Shared/OrthographicProjection.cs @@ -44,8 +44,8 @@ namespace MapControl } return new Point( - Wgs84EquatorialRadius * Math.Cos(lat) * Math.Sin(dLon), - Wgs84EquatorialRadius * (Math.Cos(lat0) * Math.Sin(lat) - Math.Sin(lat0) * Math.Cos(lat) * Math.Cos(dLon))); + Wgs84MeanRadius * Math.Cos(lat) * Math.Sin(dLon), + Wgs84MeanRadius * (Math.Cos(lat0) * Math.Sin(lat) - Math.Sin(lat0) * Math.Cos(lat) * Math.Cos(dLon))); } public override Location MapToLocation(double x, double y) @@ -55,8 +55,8 @@ namespace MapControl return new Location(Center.Latitude, Center.Longitude); } - x /= Wgs84EquatorialRadius; - y /= Wgs84EquatorialRadius; + x /= Wgs84MeanRadius; + y /= Wgs84MeanRadius; var r2 = x * x + y * y; if (r2 > 1d) diff --git a/MapControl/Shared/StereographicProjection.cs b/MapControl/Shared/StereographicProjection.cs index d4a3719b..ef3e396d 100644 --- a/MapControl/Shared/StereographicProjection.cs +++ b/MapControl/Shared/StereographicProjection.cs @@ -36,7 +36,7 @@ namespace MapControl Center.GetAzimuthDistance(latitude, longitude, out double azimuth, out double distance); - var mapDistance = Math.Tan(distance / 2d) * 2d * Wgs84EquatorialRadius; + var mapDistance = Math.Tan(distance / 2d) * 2d * Wgs84MeanRadius; return new Point(mapDistance * Math.Sin(azimuth), mapDistance * Math.Cos(azimuth)); } @@ -50,7 +50,7 @@ namespace MapControl var azimuth = Math.Atan2(x, y); var mapDistance = Math.Sqrt(x * x + y * y); - var distance = 2d * Math.Atan(mapDistance / (2d * Wgs84EquatorialRadius)); + var distance = 2d * Math.Atan(mapDistance / (2d * Wgs84MeanRadius)); return Center.GetLocation(azimuth, distance); } diff --git a/MapControl/Shared/WorldMercatorProjection.cs b/MapControl/Shared/WorldMercatorProjection.cs index 0b934c58..72abbbdf 100644 --- a/MapControl/Shared/WorldMercatorProjection.cs +++ b/MapControl/Shared/WorldMercatorProjection.cs @@ -13,6 +13,8 @@ namespace MapControl /// public class WorldMercatorProjection : MapProjection { + public static readonly double Wgs84Eccentricity = Math.Sqrt((2d - Wgs84Flattening) * Wgs84Flattening); + public const string DefaultCrsId = "EPSG:3395"; public WorldMercatorProjection() diff --git a/MapProjections/Shared/WorldMercatorProjection.cs b/MapProjections/Shared/WorldMercatorProjection.cs index 3019e4d9..874f69b0 100644 --- a/MapProjections/Shared/WorldMercatorProjection.cs +++ b/MapProjections/Shared/WorldMercatorProjection.cs @@ -43,7 +43,7 @@ namespace MapControl.Projections public override Point RelativeScale(double latitude, double longitude) { var lat = latitude * Math.PI / 180d; - var eSinLat = Wgs84Eccentricity * Math.Sin(lat); + var eSinLat = MapControl.WorldMercatorProjection.Wgs84Eccentricity * Math.Sin(lat); var k = Math.Sqrt(1d - eSinLat * eSinLat) / Math.Cos(lat); // p.44 (7-8) return new Point(k, k);