Moved EquatorialRadius property to MapProjection

This commit is contained in:
ClemensFischer 2026-01-24 14:58:16 +01:00
parent aaa69a3c03
commit d71d6c6e01
14 changed files with 53 additions and 70 deletions

View file

@ -63,13 +63,13 @@ namespace MapControl
var c = Math.Acos(p.CosC);
var k = c / Math.Sin(c); // p.195 (25-2)
return new Point(EarthRadius * k * p.X, EarthRadius * k * p.Y); // p.195 (22-4/5)
return new Point(EquatorialRadius * k * p.X, EquatorialRadius * k * p.Y); // p.195 (22-4/5)
}
public override Location MapToLocation(double x, double y)
{
var rho = Math.Sqrt(x * x + y * y);
var c = rho / EarthRadius; // p.196 (25-15)
var c = rho / EquatorialRadius; // p.196 (25-15)
return GetLocation(x, y, rho, Math.Sin(c));
}

View file

@ -5,9 +5,6 @@ using System.Windows.Media;
namespace MapControl
{
/// <summary>
/// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395), p.141.
/// </summary>
public abstract class AzimuthalProjection : MapProjection
{
protected AzimuthalProjection()
@ -16,8 +13,6 @@ namespace MapControl
Type = MapProjectionType.Azimuthal;
}
public double EarthRadius { get; set; } = Wgs84MeanRadius;
public readonly struct ProjectedPoint
{
public double X { get; }

View file

@ -35,12 +35,12 @@ namespace MapControl
public override Point? LocationToMap(double latitude, double longitude)
{
return new Point(Wgs84MeterPerDegree * longitude, Wgs84MeterPerDegree * latitude);
return new Point(MeterPerDegree * longitude, MeterPerDegree * latitude);
}
public override Location MapToLocation(double x, double y)
{
return new Location(y / Wgs84MeterPerDegree, x / Wgs84MeterPerDegree);
return new Location(y / MeterPerDegree, x / MeterPerDegree);
}
}
}

View file

@ -46,13 +46,13 @@ namespace MapControl
var k = 1d / p.CosC; // p.165 (22-3)
return new Point(EarthRadius * k * p.X, EarthRadius * k * p.Y); // p.165 (22-4/5)
return new Point(EquatorialRadius * k * p.X, EquatorialRadius * k * p.Y); // p.165 (22-4/5)
}
public override Location MapToLocation(double x, double y)
{
var rho = Math.Sqrt(x * x + y * y);
var c = Math.Atan(rho / EarthRadius); // p.167 (22-16)
var c = Math.Atan(rho / EquatorialRadius); // p.167 (22-16)
return GetLocation(x, y, rho, Math.Sin(c));
}

View file

@ -17,6 +17,11 @@ namespace MapControl
#endif
public class Location : IEquatable<Location>
{
// Arithmetic mean radius (2*a + b) / 3 == (1 - f/3) * a.
// See https://en.wikipedia.org/wiki/Earth_radius#Arithmetic_mean_radius.
//
public const double Wgs84MeanRadius = (1d - MapProjection.Wgs84Flattening / 3d) * MapProjection.Wgs84EquatorialRadius;
public Location()
{
}
@ -75,7 +80,7 @@ namespace MapControl
/// <summary>
/// Calculates great circle azimuth in degrees and distance in meters between this and the specified Location.
/// </summary>
public (double, double) GetAzimuthDistance(Location location, double earthRadius = MapProjection.Wgs84MeanRadius)
public (double, double) GetAzimuthDistance(Location location, double earthRadius = Wgs84MeanRadius)
{
var lat1 = Latitude * Math.PI / 180d;
var lon1 = Longitude * Math.PI / 180d;
@ -100,7 +105,7 @@ namespace MapControl
/// <summary>
/// Calculates great distance in meters between this and the specified Location.
/// </summary>
public double GetDistance(Location location, double earthRadius = MapProjection.Wgs84MeanRadius)
public double GetDistance(Location location, double earthRadius = Wgs84MeanRadius)
{
(var _, var distance) = GetAzimuthDistance(location, earthRadius);
@ -110,7 +115,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.Wgs84MeanRadius)
public Location GetLocation(double azimuth, double distance, double earthRadius = Wgs84MeanRadius)
{
var lat1 = Latitude * Math.PI / 180d;
var lon1 = Longitude * Math.PI / 180d;

View file

@ -31,11 +31,6 @@ namespace MapControl
public const double Wgs84Flattening = 1d / 298.257223563;
public const double Wgs84MeterPerDegree = Wgs84EquatorialRadius * Math.PI / 180d;
// Arithmetic mean radius (2*a + b) / 3 == (1 - f/3) * a.
// See https://en.wikipedia.org/wiki/Earth_radius#Arithmetic_mean_radius.
//
public const double Wgs84MeanRadius = (1d - Wgs84Flattening / 3d) * Wgs84EquatorialRadius;
public static MapProjectionFactory Factory
{
get => field ??= new MapProjectionFactory();
@ -60,6 +55,13 @@ namespace MapControl
/// </summary>
public string CrsId { get; protected set; } = "";
/// <summary>
/// The earth ellipsoid semi-major axis, or spherical earth radius, in meters.
/// </summary>
public double EquatorialRadius { get; set; } = Wgs84EquatorialRadius;
public double MeterPerDegree => EquatorialRadius * Math.PI / 180d;
private Location center;
private bool updateCenter = hasCenter;

View file

@ -42,13 +42,13 @@ namespace MapControl
return null;
}
return new Point(EarthRadius * p.X, EarthRadius * p.Y); // p.149 (20-3/4)
return new Point(EquatorialRadius * p.X, EquatorialRadius * p.Y); // p.149 (20-3/4)
}
public override Location MapToLocation(double x, double y)
{
var rho = Math.Sqrt(x * x + y * y);
var sinC = rho / EarthRadius; // p.150 (20-19)
var sinC = rho / EquatorialRadius; // p.150 (20-19)
return GetLocation(x, y, rho, sinC);
}

View file

@ -9,8 +9,8 @@ using Avalonia;
namespace MapControl
{
/// <summary>
/// Elliptical Polar Stereographic Projection with a given scale factor at the pole and
/// optional false easting and northing, as used by the UPS North and UPS South projections.
/// Elliptical Polar Stereographic Projection with scale factor at the pole and
/// false easting and northing, as used by the UPS North and UPS South projections.
/// See "Map Projections - A Working Manual" (https://pubs.usgs.gov/publication/pp1395), p.154-163.
/// </summary>
public class PolarStereographicProjection : MapProjection
@ -20,7 +20,6 @@ namespace MapControl
Type = MapProjectionType.Azimuthal;
}
public double EquatorialRadius { get; set; } = Wgs84EquatorialRadius;
public double Flattening { get; set; } = Wgs84Flattening;
public double ScaleFactor { get; set; } = 0.994;
public double FalseEasting { get; set; } = 2e6;

View file

@ -39,13 +39,13 @@ namespace MapControl
var p = GetProjectedPoint(latitude, longitude);
var k = 2d / (1d + p.CosC); // p.157 (21-4), k0 == 1
return new Point(EarthRadius * k * p.X, EarthRadius * k * p.Y); // p.157 (21-2/3)
return new Point(EquatorialRadius * k * p.X, EquatorialRadius * k * p.Y); // p.157 (21-2/3)
}
public override Location MapToLocation(double x, double y)
{
var rho = Math.Sqrt(x * x + y * y);
var c = 2d * Math.Atan(rho / (2d * EarthRadius)); // p.159 (21-15), k0 == 1
var c = 2d * Math.Atan(rho / (2d * EquatorialRadius)); // p.159 (21-15), k0 == 1
return GetLocation(x, y, rho, Math.Sin(c));
}

View file

@ -49,7 +49,6 @@ namespace MapControl
}
}
public double EquatorialRadius { get; set; } = Wgs84EquatorialRadius;
public double ScaleFactor { get; set; } = 0.9996;
public double CentralMeridian { get; set; }
public double FalseEasting { get; set; }

View file

@ -21,31 +21,12 @@ namespace MapControl
Type = MapProjectionType.TransverseCylindrical;
}
public double Flattening { get; set; } = Wgs84Flattening;
public double ScaleFactor { get; set; } = 0.9996;
public double CentralMeridian { get; set; }
public double FalseEasting { get; set; }
public double FalseNorthing { get; set; }
public double EquatorialRadius
{
get;
set
{
field = value;
M0 = MeridianDistance(LatitudeOfOrigin * Math.PI / 180d);
}
} = Wgs84EquatorialRadius;
public double Flattening
{
get;
set
{
field = value;
M0 = MeridianDistance(LatitudeOfOrigin * Math.PI / 180d);
}
} = Wgs84Flattening;
public double LatitudeOfOrigin
{
get;
@ -62,11 +43,11 @@ namespace MapControl
var e4 = e2 * e2;
var e6 = e2 * e4;
return EquatorialRadius *
((1d - e2 / 4d - 3d / 64d * e4 - 5d / 256d * e6) * phi -
(3d / 8d * e2 + 3d / 32d * e4 + 45d / 1024d * e6) * Math.Sin(2d * phi) +
(15d / 256d * e4 + 45d / 1024d * e6) * Math.Sin(4d * phi) -
35d / 3072d * e6 * Math.Sin(6d * phi)); // (3-21)
return EquatorialRadius * (
(1d - e2 / 4d - e4 * 3d / 64d - e6 * 5d / 256d) * phi -
(e2 * 3d / 8d + e4 * 3d / 32d + e6 * 45d / 1024d) * Math.Sin(2d * phi) +
(e4 * 15d / 256d + e6 * 45d / 1024d) * Math.Sin(4d * phi) -
e6 * 35d / 3072d * Math.Sin(6d * phi)); // (3-21)
}
public override Matrix RelativeScale(double latitude, double longitude)
@ -146,12 +127,12 @@ namespace MapControl
var e14 = e1 * e13;
var M = M0 + (y - FalseNorthing) / ScaleFactor; // (8-20)
var mu = M / (EquatorialRadius * (1d - e2 / 4d - 3d * e4 / 64d - 5d * e6 / 256d)); // (7-19)
var mu = M / (EquatorialRadius * (1d - e2 / 4d - e4 * 3d / 64d - e6 * 5d / 256d)); // (7-19)
var phi1 = mu +
(3d * e1 / 2d - 27d * e13 / 32d) * Math.Sin(2d * mu) +
(21d * e12 / 16d - 55d * e14 / 32d) * Math.Sin(4d * mu) +
151d * e13 / 96d * Math.Sin(6d * mu) +
1097d * e14 / 512d * Math.Sin(8d * mu); // (3-26)
(e1 * 3d / 2d - e13 * 27d / 32d) * Math.Sin(2d * mu) +
(e12 * 21d / 16d - e14 * 55d / 32d) * Math.Sin(4d * mu) +
e13 * 151d / 96d * Math.Sin(6d * mu) +
e14 * 1097d / 512d * Math.Sin(8d * mu); // (3-26)
var sinPhi1 = Math.Sin(phi1);
var cosPhi1 = Math.Cos(phi1);

View file

@ -37,15 +37,15 @@ namespace MapControl
public override Point? LocationToMap(double latitude, double longitude)
{
return new Point(
Wgs84MeterPerDegree * longitude,
Wgs84MeterPerDegree * LatitudeToY(latitude));
MeterPerDegree * longitude,
MeterPerDegree * LatitudeToY(latitude));
}
public override Location MapToLocation(double x, double y)
{
return new Location(
YToLatitude(y / Wgs84MeterPerDegree),
x / Wgs84MeterPerDegree);
YToLatitude(y / MeterPerDegree),
x / MeterPerDegree);
}
public static double LatitudeToY(double latitude)

View file

@ -39,15 +39,15 @@ namespace MapControl
public override Point? LocationToMap(double latitude, double longitude)
{
return new Point(
Wgs84MeterPerDegree * longitude,
Wgs84MeterPerDegree * LatitudeToY(latitude));
MeterPerDegree * longitude,
MeterPerDegree * LatitudeToY(latitude));
}
public override Location MapToLocation(double x, double y)
{
return new Location(
YToLatitude(y / Wgs84MeterPerDegree),
x / Wgs84MeterPerDegree);
YToLatitude(y / MeterPerDegree),
x / MeterPerDegree);
}
public static double RelativeScale(double latitude)
@ -93,10 +93,10 @@ namespace MapControl
var chi = Math.PI / 2d - 2d * Math.Atan(t); // p.45 (7-13)
return chi +
(e2 / 2d + 5d * e4 / 24d + e6 / 12d + 13d * e8 / 360d) * Math.Sin(2d * chi) +
(7d * e4 / 48d + 29d * e6 / 240d + 811d * e8 / 11520d) * Math.Sin(4d * chi) +
(7d * e6 / 120d + 81d * e8 / 1120d) * Math.Sin(6d * chi) +
4279d * e8 / 161280d * Math.Sin(8d * chi); // p.45 (3-5)
(e2 / 2d + e4 * 5d / 24d + e6 / 12d + e8 * 13d / 360d) * Math.Sin(2d * chi) +
(e4 * 7d / 48d + e6 * 29d / 240d + e8 * 811d / 11520d) * Math.Sin(4d * chi) +
(e6 * 7d / 120d + e8 * 81d / 1120d) * Math.Sin(6d * chi) +
e8 * 4279d / 161280d * Math.Sin(8d * chi); // p.45 (3-5)
}
}
}

View file

@ -211,8 +211,10 @@
Description="© [BKG](https://gdz.bkg.bund.de/index.php/default/webdienste/topplus-produkte/wms-topplusopen-mit-layer-fur-normalausgabe-und-druck-wms-topplus-open.html)"/>
</tools:MapLayerMenuItem>
<tools:MapLayerMenuItem Text="OpenStreetMap WMS">
<!--http://t400:8080/geoserver/terrestris/wms?service=WMS&version=1.1.0&request=GetMap&layers=terrestris%3AOSM-Overlay-WMS&bbox=-2.0037508342789244E7%2C-2.5819498513543323E7%2C2.0037508342789244E7%2C2.58194985135433E7&width=596&height=768&srs=EPSG%3A3857&styles=&format=application/openlayers-->
<map:WmsImageLayer
ServiceUri="https://ows.terrestris.de/osm/service"
ServiceUri="http://t400:8080/geoserver/terrestris/wms"
RequestLayers="OSM-WMS"
Description="© [terrestris GmbH &amp; Co. KG](http://ows.terrestris.de/) © [OpenStreetMap contributors](http://www.openstreetmap.org/copyright)"/>
</tools:MapLayerMenuItem>