Updated StereographicProjection

This commit is contained in:
ClemensFischer 2026-02-01 10:53:12 +01:00
parent cfed3575ac
commit a4bd11e26d
2 changed files with 40 additions and 27 deletions

View file

@ -27,10 +27,25 @@ namespace MapControl
public override double GridConvergence(double latitude, double longitude)
{
var p0 = LocationToMap(latitude, longitude);
var p1 = LocationToMap(latitude + 1e-3, longitude);
var phi0 = LatitudeOfOrigin * Math.PI / 180d; // φ1
var phi1 = latitude * Math.PI / 180d;
var phi2 = (latitude + 1e-3) * Math.PI / 180d;
var dLambda = (longitude - CentralMeridian) * Math.PI / 180d; // λ - λ0
var sinPhi0 = Math.Sin(phi0);
var cosPhi0 = Math.Cos(phi0);
var sinPhi1 = Math.Sin(phi1);
var cosPhi1 = Math.Cos(phi1);
var sinPhi2 = Math.Sin(phi2);
var cosPhi2 = Math.Cos(phi2);
var sinLambda = Math.Sin(dLambda);
var cosLambda = Math.Cos(dLambda);
var k1 = 2d / (1d + sinPhi0 * sinPhi1 + cosPhi0 * cosPhi1 * cosLambda);
var k2 = 2d / (1d + sinPhi0 * sinPhi2 + cosPhi0 * cosPhi2 * cosLambda);
var dCosPhi = k2 * cosPhi2 - k1 * cosPhi1;
var dSinPhi = k2 * sinPhi2 - k1 * sinPhi1;
return Math.Atan2(p0.X - p1.X, p1.Y - p0.Y) * 180d / Math.PI;
return Math.Atan2(-sinLambda * dCosPhi,
cosPhi0 * dSinPhi - sinPhi0 * cosLambda * dCosPhi) * 180d / Math.PI;
}
public override Matrix RelativeTransform(double latitude, double longitude)
@ -40,20 +55,18 @@ namespace MapControl
public override Point LocationToMap(double latitude, double longitude)
{
var phi = latitude * Math.PI / 180d;
var phi1 = LatitudeOfOrigin * Math.PI / 180d;
var phi0 = LatitudeOfOrigin * Math.PI / 180d; // φ1
var phi = latitude * Math.PI / 180d; // φ
var dLambda = (longitude - CentralMeridian) * Math.PI / 180d; // λ - λ0
var cosPhi = Math.Cos(phi);
var sinPhi0 = Math.Sin(phi0);
var cosPhi0 = Math.Cos(phi0);
var sinPhi = Math.Sin(phi);
var cosPhi1 = Math.Cos(phi1);
var sinPhi1 = Math.Sin(phi1);
var cosLambda = Math.Cos(dLambda);
var cosPhi = Math.Cos(phi);
var sinLambda = Math.Sin(dLambda);
var cosPhiCosLambda = cosPhi * Math.Cos(dLambda);
var x = cosPhi * sinLambda;
var y = cosPhi1 * sinPhi - sinPhi1 * cosPhi * cosLambda;
var cosC = sinPhi1 * sinPhi + cosPhi1 * cosPhi * cosLambda; // (5-3)
cosC = Math.Min(Math.Max(cosC, -1d), 1d); // protect against rounding errors
var k = 2d / (1d + cosC); // p.157 (21-4), k0 == 1
var y = cosPhi0 * sinPhi - sinPhi0 * cosPhiCosLambda;
var k = 2d / (1d + sinPhi0 * sinPhi + cosPhi0 * cosPhiCosLambda); // p.157 (21-4), k0 == 1
return new Point(
EquatorialRadius * k * x,
@ -67,10 +80,10 @@ namespace MapControl
var cosC = Math.Cos(c);
var sinC = Math.Sin(c);
var phi1 = LatitudeOfOrigin * Math.PI / 180d;
var cosPhi1 = Math.Cos(phi1);
var sinPhi1 = Math.Sin(phi1);
var phi = Math.Asin(cosC * sinPhi1 + y * sinC * cosPhi1 / rho); // (20-14)
var phi0 = LatitudeOfOrigin * Math.PI / 180d; // φ1
var cosPhi0 = Math.Cos(phi0);
var sinPhi0 = Math.Sin(phi0);
var phi = Math.Asin(cosC * sinPhi0 + y * sinC * cosPhi0 / rho); // (20-14)
double u, v;
if (LatitudeOfOrigin == 90d) // (20-16)
@ -86,7 +99,7 @@ namespace MapControl
else // (20-15)
{
u = x * sinC;
v = rho * cosPhi1 * cosC - y * sinPhi1 * sinC;
v = rho * cosPhi0 * cosC - y * sinPhi0 * sinC;
}
return new Location(