Map projections

This commit is contained in:
ClemensFischer 2025-02-12 19:46:25 +01:00
parent e90bdcf371
commit 7c8393d785
7 changed files with 24 additions and 27 deletions

View file

@ -32,7 +32,7 @@ namespace MapControl
EquatorialRadius = 6378137d;
Flattening = 1d / 298.257222101;
ScaleFactor = 0.9996;
CentralMeridian = Zone * 6d - 183d;
CentralMeridian = zone * 6d - 183d;
FalseEasting = 5e5;
FalseNorthing = 0d;
}

View file

@ -32,7 +32,7 @@ namespace MapControl
EquatorialRadius = 6378206.4;
Flattening = 1d / 294.978698213898;
ScaleFactor = 0.9996;
CentralMeridian = Zone * 6d - 183d;
CentralMeridian = zone * 6d - 183d;
FalseEasting = 5e5;
FalseNorthing = 0d;
}

View file

@ -32,7 +32,7 @@ namespace MapControl
EquatorialRadius = 6378137d;
Flattening = 1d / 298.257222101;
ScaleFactor = 0.9996;
CentralMeridian = Zone * 6d - 183d;
CentralMeridian = zone * 6d - 183d;
FalseEasting = 5e5;
FalseNorthing = 0d;
}

View file

@ -94,12 +94,11 @@ namespace MapControl
}
var e = Math.Sqrt((2d - Flattening) * Flattening);
var r = Math.Sqrt(x * x + y * y); // p.162 (20-18)
var t = r * Math.Sqrt(Math.Pow(1d + e, 1d + e) * Math.Pow(1d - e, 1d - e))
/ (2d * EquatorialRadius * ScaleFactor); // p.162 (21-39)
var lat = WorldMercatorProjection.LatitudeFromSeriesApproximation(e, t); // p.162 (3-5)
var lat = WorldMercatorProjection.ApproximateLatitude(e, t); // p.162 (3-5)
var lon = Math.Atan2(x, -y); // p.162 (20-16)
if (!IsNorth)

View file

@ -15,12 +15,6 @@ namespace MapControl
/// </summary>
public class TransverseMercatorProjection : MapProjection
{
#if NETFRAMEWORK || UWP
private static double Math_Atanh(double x) => Math.Log((1d + x) / (1d - x)) / 2d;
#else
private static double Math_Atanh(double x) => Math.Atanh(x);
#endif
public double EquatorialRadius { get; set; } = Wgs84EquatorialRadius;
public double Flattening { get; set; } = Wgs84Flattening;
public double ScaleFactor { get; set; } = 0.9996;
@ -35,19 +29,22 @@ namespace MapControl
public override Point GetRelativeScale(Location location)
{
// Precise enough, avoid calculations as in LocationToMap.
//
return new Point(ScaleFactor, ScaleFactor);
}
public override Point? LocationToMap(Location location)
{
#if NETFRAMEWORK || UWP
double Atanh(double x) => Math.Log((1d + x) / (1d - x)) / 2d;
#else
static double Atanh(double x) => Math.Atanh(x);
#endif
var n = Flattening / (2d - Flattening);
var n2 = n * n;
var n3 = n * n2;
var k0A = ScaleFactor * EquatorialRadius / (1d + n) * (1d + n2 / 4d + n2 * n2 / 64d);
// α_j
var alpha1 = n / 2d - n2 * 2d / 3d + n3 * 5d / 16d;
var alpha2 = n2 * 13d / 48d - n3 * 3d / 5d;
var alpha3 = n3 * 61d / 240d;
@ -59,15 +56,14 @@ namespace MapControl
var lambda = (location.Longitude - CentralMeridian) * Math.PI / 180d;
var s = 2d * Math.Sqrt(n) / (1d + n);
var sinLat = Math.Sin(phi);
var t = Math.Sinh(Math_Atanh(sinLat) - s * Math_Atanh(s * sinLat));
var sinPhi = Math.Sin(phi);
var t = Math.Sinh(Atanh(sinPhi) - s * Atanh(s * sinPhi));
// ξ'
var xi_ = Math.Atan(t / Math.Cos(lambda));
// η'
var eta_ = Math_Atanh(Math.Sin(lambda) / Math.Sqrt(1d + t * t));
var eta_ = Atanh(Math.Sin(lambda) / Math.Sqrt(1d + t * t));
// ξ
var xi = xi_
@ -91,13 +87,14 @@ namespace MapControl
var n = Flattening / (2d - Flattening);
var n2 = n * n;
var n3 = n * n2;
var k0A = ScaleFactor * EquatorialRadius / (1d + n) * (1d + n2 / 4d + n2 * n2 / 64d);
// β_j
var beta1 = n / 2d - n2 * 2d / 3d + n3 * 37d / 96d;
var beta2 = n2 / 48d + n3 / 15d;
var beta3 = n3 * 17d / 480d;
// δ_j
var delta1 = n * 2d - n2 * 2d / 3d - n3 * 2d;
var delta2 = n2 * 7d / 3d - n3 * 8d / 5d;
var delta3 = n3 * 56d / 15d;

View file

@ -24,6 +24,11 @@ namespace MapControl
public Wgs84UtmProjection(int zone, bool north)
{
SetZone(zone, north);
EquatorialRadius = Wgs84EquatorialRadius;
Flattening = Wgs84Flattening;
ScaleFactor = 0.9996;
FalseEasting = 5e5;
}
protected void SetZone(int zone, bool north)
@ -36,12 +41,8 @@ namespace MapControl
Zone = zone;
IsNorth = north;
CrsId = $"EPSG:{(north ? 32600 : 32700) + zone}";
EquatorialRadius = Wgs84EquatorialRadius;
Flattening = Wgs84Flattening;
ScaleFactor = 0.9996;
CentralMeridian = Zone * 6d - 183d;
FalseEasting = 5e5;
FalseNorthing = IsNorth ? 0d : 1e7;
CentralMeridian = zone * 6d - 183d;
FalseNorthing = north ? 0d : 1e7;
}
}

View file

@ -75,10 +75,10 @@ namespace MapControl
{
var t = Math.Exp(-y * Math.PI / 180d); // p.44 (7-10)
return LatitudeFromSeriesApproximation(Wgs84Eccentricity, t) * 180d / Math.PI;
return ApproximateLatitude(Wgs84Eccentricity, t) * 180d / Math.PI;
}
internal static double LatitudeFromSeriesApproximation(double e, double t)
internal static double ApproximateLatitude(double e, double t)
{
var e_2 = e * e;
var e_4 = e_2 * e_2;