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);