Use mean radius in distance calculations

This commit is contained in:
ClemensFischer 2026-01-09 08:13:07 +01:00
parent c187f323a8
commit 23a8e49efb
8 changed files with 20 additions and 14 deletions

View file

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

View file

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

View file

@ -115,7 +115,7 @@ namespace MapControl
/// <summary>
/// Calculates the great circle distance in meters between this and the specified Location.
/// </summary>
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
/// <summary>
/// Calculates the Location on a great circle at the specified azimuth in degrees and distance in meters from this Location.
/// </summary>
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);
}

View file

@ -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
{

View file

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

View file

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

View file

@ -13,6 +13,8 @@ namespace MapControl
/// </summary>
public class WorldMercatorProjection : MapProjection
{
public static readonly double Wgs84Eccentricity = Math.Sqrt((2d - Wgs84Flattening) * Wgs84Flattening);
public const string DefaultCrsId = "EPSG:3395";
public WorldMercatorProjection()

View file

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