Updated TransverseMercatorProjection projections

This commit is contained in:
ClemensFischer 2026-01-23 20:37:17 +01:00
parent e475653a76
commit 3d000ce54b
2 changed files with 57 additions and 23 deletions

View file

@ -63,13 +63,6 @@ namespace MapControl
public override Point? LocationToMap(double latitude, double longitude)
{
#if NETFRAMEWORK
static double Atanh(double x) => Math.Log((1d + x) / (1d - x)) / 2d;
#else
static double Atanh(double x) => Math.Atanh(x);
#endif
// k0 * A
var k0A = ScaleFactor * EquatorialRadius * f1;
// φ
var phi = latitude * Math.PI / 180d;
var sinPhi = Math.Sin(phi);
@ -81,20 +74,20 @@ namespace MapControl
var xi_ = Math.Atan2(t, Math.Cos(dLambda));
// η'
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 eta = eta_ +
// k0 * A
var k0A = ScaleFactor * EquatorialRadius * f1;
// E
var x = FalseEasting + k0A * (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_);
a3 * Math.Cos(6d * xi_) * Math.Sinh(6d * eta_));
// N
var y = FalseNorthing + k0A * (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_));
return new Point(
k0A * eta + FalseEasting,
k0A * xi + FalseNorthing);
return new Point(x, y);
}
public override Location MapToLocation(double x, double y)
@ -129,5 +122,46 @@ namespace MapControl
phi * 180d / Math.PI,
dLambda * 180d / Math.PI + CentralMeridian);
}
#if NETFRAMEWORK
private static double Atanh(double x) => Math.Log((1d + x) / (1d - x)) / 2d;
#else
private static double Atanh(double x) => Math.Atanh(x);
#endif
/*
* Relative scale is usually < 1.001 and hence neglectable.
*
public override Matrix RelativeScale(double latitude, double longitude)
{
// φ
var phi = latitude * Math.PI / 180d;
var sinPhi = Math.Sin(phi);
// t
var t = Math.Sinh(Atanh(sinPhi) - f2 * Atanh(f2 * sinPhi));
// λ - λ0
var dLambda = (longitude - CentralMeridian) * Math.PI / 180d;
var cosLambda = Math.Cos(dLambda);
// ξ'
var xi_ = Math.Atan2(t, cosLambda);
// η'
var eta_ = Atanh(Math.Sin(dLambda) / Math.Sqrt(1d + t * t));
// σ
var sigma = 1 +
2d * a1 * Math.Cos(2d * xi_) * Math.Cosh(2d * eta_) +
4d * a2 * Math.Cos(4d * xi_) * Math.Cosh(4d * eta_) +
6d * a3 * Math.Cos(6d * xi_) * Math.Cosh(6d * eta_);
// τ
var tau =
2d * a1 * Math.Sin(2d * xi_) * Math.Sinh(2d * eta_) +
4d * a2 * Math.Sin(4d * xi_) * Math.Sinh(4d * eta_) +
6d * a3 * Math.Sin(6d * xi_) * Math.Sinh(6d * eta_);
// h = k/k0 for relative scale
var n = Flattening / (2d - Flattening);
var u = (1d - n) / (1d + n) * Math.Tan(phi);
var h = f1 * Math.Sqrt((1d + u * u) * (sigma * sigma + tau * tau) / (t * t + cosLambda * cosLambda));
return new Matrix(h, 0d, 0d, h, 0d, 0d);
}*/
}
}

View file

@ -63,15 +63,15 @@ namespace MapControl
var e6 = e2 * e4;
return EquatorialRadius *
((1d - e2 / 4d - 3d * e4 / 64d - 5d * e6 / 256d) * phi -
(3d * e2 / 8d + 3d * e4 / 32d + 45d * e6 / 1024d) * Math.Sin(2d * phi) +
(15d * e4 / 256d + 45d * e6 / 1024d) * Math.Sin(4d * phi) -
35d * e6 / 3072d * Math.Sin(6d * phi)); // (3-21)
((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)
}
public override Matrix RelativeScale(double latitude, double longitude)
{
var k = ScaleFactor;
var k = 1d; // omit k0 for relative scale
if (latitude > -90d && latitude < 90d)
{