diff --git a/MapControl/Shared/Etrs89UtmProjection.cs b/MapControl/Shared/Etrs89UtmProjection.cs index 9c9d6a5d..fbd6316d 100644 --- a/MapControl/Shared/Etrs89UtmProjection.cs +++ b/MapControl/Shared/Etrs89UtmProjection.cs @@ -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; } diff --git a/MapControl/Shared/Nad27UtmProjection.cs b/MapControl/Shared/Nad27UtmProjection.cs index 32f707c9..27dc4c4d 100644 --- a/MapControl/Shared/Nad27UtmProjection.cs +++ b/MapControl/Shared/Nad27UtmProjection.cs @@ -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; } diff --git a/MapControl/Shared/Nad83UtmProjection.cs b/MapControl/Shared/Nad83UtmProjection.cs index d73c5e39..a6343294 100644 --- a/MapControl/Shared/Nad83UtmProjection.cs +++ b/MapControl/Shared/Nad83UtmProjection.cs @@ -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; } diff --git a/MapControl/Shared/PolarStereographicProjection.cs b/MapControl/Shared/PolarStereographicProjection.cs index bc28dfcc..43107a4b 100644 --- a/MapControl/Shared/PolarStereographicProjection.cs +++ b/MapControl/Shared/PolarStereographicProjection.cs @@ -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) diff --git a/MapControl/Shared/TransverseMercatorProjection.cs b/MapControl/Shared/TransverseMercatorProjection.cs index e34f557d..3057094d 100644 --- a/MapControl/Shared/TransverseMercatorProjection.cs +++ b/MapControl/Shared/TransverseMercatorProjection.cs @@ -15,12 +15,6 @@ namespace MapControl /// 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; diff --git a/MapControl/Shared/Wgs84UtmProjection.cs b/MapControl/Shared/Wgs84UtmProjection.cs index 0697f1bc..77e7e216 100644 --- a/MapControl/Shared/Wgs84UtmProjection.cs +++ b/MapControl/Shared/Wgs84UtmProjection.cs @@ -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; } } diff --git a/MapControl/Shared/WorldMercatorProjection.cs b/MapControl/Shared/WorldMercatorProjection.cs index 62f1dd37..f74d5a7d 100644 --- a/MapControl/Shared/WorldMercatorProjection.cs +++ b/MapControl/Shared/WorldMercatorProjection.cs @@ -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;