Added AzimuthalProjection.EarthRadius

This commit is contained in:
ClemensFischer 2026-01-16 20:22:45 +01:00
parent 60152fef93
commit d825b50062
9 changed files with 63 additions and 60 deletions

View file

@ -34,7 +34,7 @@ namespace MapControl
Center.GetAzimuthDistance(latitude, longitude, out double azimuth, out double distance);
var mapDistance = distance * Wgs84MeanRadius;
var mapDistance = distance * EarthRadius;
return new Point(mapDistance * Math.Sin(azimuth), mapDistance * Math.Cos(azimuth));
}
@ -48,7 +48,7 @@ namespace MapControl
var azimuth = Math.Atan2(x, y);
var mapDistance = Math.Sqrt(x * x + y * y);
var distance = mapDistance / Wgs84MeanRadius;
var distance = mapDistance / EarthRadius;
return Center.GetLocation(azimuth, distance);
}

View file

@ -7,7 +7,7 @@ using Avalonia;
namespace MapControl
{
/// <summary>
/// Base class for azimuthal map projections.
/// Base class for spherical azimuthal map projections.
/// </summary>
public abstract class AzimuthalProjection : MapProjection
{
@ -16,6 +16,8 @@ namespace MapControl
Type = MapProjectionType.Azimuthal;
}
public double EarthRadius { get; set; } = Wgs84MeanRadius;
public override Rect? BoundingBoxToMap(BoundingBox boundingBox)
{
Rect? rect = null;

View file

@ -35,7 +35,7 @@ namespace MapControl
Center.GetAzimuthDistance(latitude, longitude, out double azimuth, out double distance);
var mapDistance = distance < Math.PI / 2d
? Math.Tan(distance) * Wgs84MeanRadius
? Math.Tan(distance) * EarthRadius
: double.PositiveInfinity;
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 = Math.Atan(mapDistance / Wgs84MeanRadius);
var distance = Math.Atan(mapDistance / EarthRadius);
return Center.GetLocation(azimuth, distance);
}

View file

@ -41,9 +41,9 @@ namespace MapControl
return null;
}
var x = Wgs84MeanRadius * Math.Cos(phi) * Math.Sin(dLambda); // p.149 (20-3)
var y = Wgs84MeanRadius * (Math.Cos(phi1) * Math.Sin(phi) -
Math.Sin(phi1) * Math.Cos(phi) * Math.Cos(dLambda)); // p.149 (20-4)
var x = EarthRadius * Math.Cos(phi) * Math.Sin(dLambda); // p.149 (20-3)
var y = EarthRadius * (Math.Cos(phi1) * Math.Sin(phi) -
Math.Sin(phi1) * Math.Cos(phi) * Math.Cos(dLambda)); // p.149 (20-4)
return new Point(x, y);
}
@ -54,8 +54,8 @@ namespace MapControl
return new Location(Center.Latitude, Center.Longitude);
}
x /= Wgs84MeanRadius;
y /= Wgs84MeanRadius;
x /= EarthRadius;
y /= EarthRadius;
var r2 = x * x + y * y;
if (r2 > 1d)

View file

@ -34,7 +34,7 @@ namespace MapControl
Center.GetAzimuthDistance(latitude, longitude, out double azimuth, out double distance);
var mapDistance = Math.Tan(distance / 2d) * 2d * Wgs84MeanRadius;
var mapDistance = Math.Tan(distance / 2d) * 2d * EarthRadius;
return new Point(mapDistance * Math.Sin(azimuth), mapDistance * Math.Cos(azimuth));
}
@ -48,7 +48,7 @@ namespace MapControl
var azimuth = Math.Atan2(x, y);
var mapDistance = Math.Sqrt(x * x + y * y);
var distance = 2d * Math.Atan(mapDistance / (2d * Wgs84MeanRadius));
var distance = 2d * Math.Atan(mapDistance / (2d * EarthRadius));
return Center.GetLocation(azimuth, distance);
}

View file

@ -87,15 +87,15 @@ namespace MapControl
// η'
var eta_ = Atanh(Math.Sin(dLambda) / Math.Sqrt(1d + t * t));
// ξ
var xi = xi_
+ a1 * Math.Sin(2d * xi_) * Math.Cosh(2d * eta_)
+ a2 * Math.Sin(4d * xi_) * Math.Cosh(4d * eta_)
+ a3 * Math.Sin(6d * xi_) * Math.Cosh(6d * eta_);
var xi = xi_ +
a1 * Math.Sin(2d * xi_) * Math.Cosh(2d * eta_) +
a2 * Math.Sin(4d * xi_) * Math.Cosh(4d * eta_) +
a3 * Math.Sin(6d * xi_) * Math.Cosh(6d * eta_);
// η
var eta = eta_
+ a1 * Math.Cos(2d * xi_) * Math.Sinh(2d * eta_)
+ a2 * Math.Cos(4d * xi_) * Math.Sinh(4d * eta_)
+ a3 * Math.Cos(6d * xi_) * Math.Sinh(6d * eta_);
var eta = eta_ +
a1 * Math.Cos(2d * xi_) * Math.Sinh(2d * eta_) +
a2 * Math.Cos(4d * xi_) * Math.Sinh(4d * eta_) +
a3 * Math.Cos(6d * xi_) * Math.Sinh(6d * eta_);
return new Point(
k0A * eta + FalseEasting,
@ -111,22 +111,22 @@ namespace MapControl
// η
var eta = (x - FalseEasting) / k0A;
// ξ'
var xi_ = xi
- b1 * Math.Sin(2d * xi) * Math.Cosh(2d * eta)
- b2 * Math.Sin(4d * xi) * Math.Cosh(4d * eta)
- b3 * Math.Sin(6d * xi) * Math.Cosh(6d * eta);
var xi_ = xi -
b1 * Math.Sin(2d * xi) * Math.Cosh(2d * eta) -
b2 * Math.Sin(4d * xi) * Math.Cosh(4d * eta) -
b3 * Math.Sin(6d * xi) * Math.Cosh(6d * eta);
// η'
var eta_ = eta
- b1 * Math.Cos(2d * xi) * Math.Sinh(2d * eta)
- b2 * Math.Cos(4d * xi) * Math.Sinh(4d * eta)
- b3 * Math.Cos(6d * xi) * Math.Sinh(6d * eta);
var eta_ = eta -
b1 * Math.Cos(2d * xi) * Math.Sinh(2d * eta) -
b2 * Math.Cos(4d * xi) * Math.Sinh(4d * eta) -
b3 * Math.Cos(6d * xi) * Math.Sinh(6d * eta);
// χ
var chi = Math.Asin(Math.Sin(xi_) / Math.Cosh(eta_));
// φ
var phi = chi
+ d1 * Math.Sin(2d * chi)
+ d2 * Math.Sin(4d * chi)
+ d3 * Math.Sin(6d * chi);
var phi = chi +
d1 * Math.Sin(2d * chi) +
d2 * Math.Sin(4d * chi) +
d3 * Math.Sin(6d * chi);
// λ - λ0
var dLambda = Math.Atan(Math.Sinh(eta_) / Math.Cos(xi_));

View file

@ -170,16 +170,17 @@ namespace MapControl
format = formatPng;
}
uriTemplate += "?Service=WMTS"
+ "&Request=GetTile"
+ "&Version=1.0.0"
+ "&Layer=" + layer
+ "&Style=" + style
+ "&Format=" + format
+ "&TileMatrixSet={TileMatrixSet}"
+ "&TileMatrix={TileMatrix}"
+ "&TileRow={TileRow}"
+ "&TileCol={TileCol}";
uriTemplate +=
"?Service=WMTS" +
"&Request=GetTile" +
"&Version=1.0.0" +
"&Layer=" + layer +
"&Style=" + style +
"&Format=" + format +
"&TileMatrixSet={TileMatrixSet}" +
"&TileMatrix={TileMatrix}" +
"&TileRow={TileRow}" +
"&TileCol={TileCol}";
}
}

View file

@ -91,11 +91,11 @@ namespace MapControl
var e8 = e2 * e6;
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)
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)
}
}
}

View file

@ -14,19 +14,19 @@ namespace MapControl.Projections
{
public WorldMercatorProjection()
{
CoordinateSystemWkt
= "PROJCS[\"WGS 84 / World Mercator\","
+ WktConstants.GeogCsWgs84 + ","
+ "PROJECTION[\"Mercator_1SP\"],"
+ "PARAMETER[\"latitude_of_origin\",0],"
+ "PARAMETER[\"central_meridian\",0],"
+ "PARAMETER[\"scale_factor\",1],"
+ "PARAMETER[\"false_easting\",0],"
+ "PARAMETER[\"false_northing\",0],"
+ "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],"
+ "AXIS[\"Easting\",EAST],"
+ "AXIS[\"Northing\",NORTH],"
+ "AUTHORITY[\"EPSG\",\"3395\"]]";
CoordinateSystemWkt =
"PROJCS[\"WGS 84 / World Mercator\"," +
WktConstants.GeogCsWgs84 + "," +
"PROJECTION[\"Mercator_1SP\"]," +
"PARAMETER[\"latitude_of_origin\",0]," +
"PARAMETER[\"central_meridian\",0]," +
"PARAMETER[\"scale_factor\",1]," +
"PARAMETER[\"false_easting\",0]," +
"PARAMETER[\"false_northing\",0]," +
"UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," +
"AXIS[\"Easting\",EAST]," +
"AXIS[\"Northing\",NORTH]," +
"AUTHORITY[\"EPSG\",\"3395\"]]";
}
public override Point RelativeScale(double latitude, double longitude)