diff --git a/MapControl/Shared/OrthographicProjection.cs b/MapControl/Shared/OrthographicProjection.cs index b4a1b7ea..85b95393 100644 --- a/MapControl/Shared/OrthographicProjection.cs +++ b/MapControl/Shared/OrthographicProjection.cs @@ -32,18 +32,19 @@ namespace MapControl return new Point(); } - var lat0 = Center.Latitude * Math.PI / 180d; - var lat = latitude * Math.PI / 180d; - var dLon = (longitude - Center.Longitude) * Math.PI / 180d; + var phi = latitude * Math.PI / 180d; + var phi1 = Center.Latitude * Math.PI / 180d; + var lambda = (longitude - Center.Longitude) * Math.PI / 180d; // λ - λ0 - if (Math.Abs(lat - lat0) > Math.PI / 2d || Math.Abs(dLon) > Math.PI / 2d) + if (Math.Abs(phi - phi1) > Math.PI / 2d || Math.Abs(lambda) > Math.PI / 2d) { return null; } - return new Point( - Wgs84MeanRadius * Math.Cos(lat) * Math.Sin(dLon), - Wgs84MeanRadius * (Math.Cos(lat0) * Math.Sin(lat) - Math.Sin(lat0) * Math.Cos(lat) * Math.Cos(dLon))); + var x = Wgs84MeanRadius * Math.Cos(phi) * Math.Sin(lambda); // p.149 (20-3) + var y = Wgs84MeanRadius * (Math.Cos(phi1) * Math.Sin(phi) - + Math.Sin(phi1) * Math.Cos(phi) * Math.Cos(lambda)); // p.149 (20-4) + return new Point(x, y); } public override Location MapToLocation(double x, double y) @@ -66,13 +67,14 @@ namespace MapControl var sinC = r; var cosC = Math.Sqrt(1 - r2); - var lat0 = Center.Latitude * Math.PI / 180d; - var cosLat0 = Math.Cos(lat0); - var sinLat0 = Math.Sin(lat0); + var phi1 = Center.Latitude * Math.PI / 180d; + var cosPhi1 = Math.Cos(phi1); + var sinPhi1 = Math.Sin(phi1); - return new Location( - 180d / Math.PI * Math.Asin(cosC * sinLat0 + y * sinC * cosLat0 / r), - 180d / Math.PI * Math.Atan2(x * sinC, r * cosC * cosLat0 - y * sinC * sinLat0) + Center.Longitude); + var phi = Math.Asin(cosC * sinPhi1 + y * sinC * cosPhi1 / r); // p.150 (20-14) + var lambda = Math.Atan2(x * sinC, r * cosC * cosPhi1 - y * sinC * sinPhi1); // p.150 (20-15) + + return new Location(180d / Math.PI * phi, 180d / Math.PI * lambda + Center.Longitude); } } } diff --git a/MapControl/Shared/PolarStereographicProjection.cs b/MapControl/Shared/PolarStereographicProjection.cs index c8c48db5..d0c01c8c 100644 --- a/MapControl/Shared/PolarStereographicProjection.cs +++ b/MapControl/Shared/PolarStereographicProjection.cs @@ -96,16 +96,16 @@ namespace MapControl 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.ApproximateLatitude(e, t); // p.162 (3-5) - var lon = Math.Atan2(x, -y); // p.162 (20-16) + var phi = WorldMercatorProjection.ApproximateLatitude(e, t); // p.162 (3-5) + var lambda = Math.Atan2(x, -y); // p.162 (20-16) if (Hemisphere == Hemisphere.South) { - lat = -lat; - lon = -lon; + phi = -phi; + lambda = -lambda; } - return new Location(lat * 180d / Math.PI, lon * 180d / Math.PI); + return new Location(phi * 180d / Math.PI, lambda * 180d / Math.PI); } } diff --git a/MapControl/Shared/WorldMercatorProjection.cs b/MapControl/Shared/WorldMercatorProjection.cs index 6ecbda1a..4f0c245d 100644 --- a/MapControl/Shared/WorldMercatorProjection.cs +++ b/MapControl/Shared/WorldMercatorProjection.cs @@ -30,9 +30,7 @@ namespace MapControl public override Point RelativeScale(double latitude, double longitude) { - var lat = latitude * Math.PI / 180d; - var eSinLat = Wgs84Eccentricity * Math.Sin(lat); - var k = Math.Sqrt(1d - eSinLat * eSinLat) / Math.Cos(lat); // p.44 (7-8) + var k = ScaleFactor(latitude); return new Point(k, k); } @@ -51,6 +49,14 @@ namespace MapControl x / Wgs84MeterPerDegree); } + public static double ScaleFactor(double latitude) + { + var phi = latitude * Math.PI / 180d; + var eSinPhi = Wgs84Eccentricity * Math.Sin(phi); + + return Math.Sqrt(1d - eSinPhi * eSinPhi) / Math.Cos(phi); // p.44 (7-8) + } + public static double LatitudeToY(double latitude) { if (latitude <= -90d) @@ -63,11 +69,11 @@ namespace MapControl return double.PositiveInfinity; } - var lat = latitude * Math.PI / 180d; - var eSinLat = Wgs84Eccentricity * Math.Sin(lat); - var f = Math.Pow((1d - eSinLat) / (1d + eSinLat), Wgs84Eccentricity / 2d); + var phi = latitude * Math.PI / 180d; + var eSinPhi = Wgs84Eccentricity * Math.Sin(phi); + var f = Math.Pow((1d - eSinPhi) / (1d + eSinPhi), Wgs84Eccentricity / 2d); - return Math.Log(Math.Tan(lat / 2d + Math.PI / 4d) * f) * 180d / Math.PI; // p.44 (7-7) + return Math.Log(Math.Tan(phi / 2d + Math.PI / 4d) * f) * 180d / Math.PI; // p.44 (7-7) } public static double YToLatitude(double y) @@ -79,18 +85,17 @@ namespace MapControl internal static double ApproximateLatitude(double e, double t) { - var e_2 = e * e; - var e_4 = e_2 * e_2; - var e_6 = e_2 * e_4; - var e_8 = e_2 * e_6; + var e2 = e * e; + var e4 = e2 * e2; + var e6 = e2 * e4; + var e8 = e2 * e6; + var chi = Math.PI / 2d - 2d * Math.Atan(t); // p.45 (7-13) - var lat = Math.PI / 2d - 2d * Math.Atan(t); // p.45 (7-13) - - return lat - + (e_2 / 2d + 5d * e_4 / 24d + e_6 / 12d + 13d * e_8 / 360d) * Math.Sin(2d * lat) - + (7d * e_4 / 48d + 29d * e_6 / 240d + 811d * e_8 / 11520d) * Math.Sin(4d * lat) - + (7d * e_6 / 120d + 81d * e_8 / 1120d) * Math.Sin(6d * lat) - + 4279d * e_8 / 161280d * Math.Sin(8d * lat); // 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) } } } diff --git a/MapProjections/Shared/GeoApiProjectionFactory.cs b/MapProjections/Shared/GeoApiProjectionFactory.cs index a9ce9c26..d7e25ec1 100644 --- a/MapProjections/Shared/GeoApiProjectionFactory.cs +++ b/MapProjections/Shared/GeoApiProjectionFactory.cs @@ -241,7 +241,7 @@ namespace MapControl.Projections { MapControl.WebMercatorProjection.DefaultCrsId => new WebMercatorProjection(), MapControl.WorldMercatorProjection.DefaultCrsId => new WorldMercatorProjection(), - MapControl.Wgs84AutoUtmProjection.DefaultCrsId => new Wgs84AutoUtmProjection(), + Wgs84AutoUtmProjection.DefaultCrsId => new Wgs84AutoUtmProjection(), _ => base.GetProjection(crsId) }; diff --git a/MapProjections/Shared/Wgs84AutoUtmProjection.cs b/MapProjections/Shared/Wgs84AutoUtmProjection.cs index 18c03ce1..e376f123 100644 --- a/MapProjections/Shared/Wgs84AutoUtmProjection.cs +++ b/MapProjections/Shared/Wgs84AutoUtmProjection.cs @@ -9,10 +9,12 @@ namespace MapControl.Projections /// public class Wgs84AutoUtmProjection : Wgs84UtmProjection { + public const string DefaultCrsId = "AUTO2:42001"; + private readonly string autoCrsId; public Wgs84AutoUtmProjection() // parameterless constructor for XAML - : this(MapControl.Wgs84AutoUtmProjection.DefaultCrsId) + : this(DefaultCrsId) { } diff --git a/MapProjections/Shared/WorldMercatorProjection.cs b/MapProjections/Shared/WorldMercatorProjection.cs index adefbc46..089e7759 100644 --- a/MapProjections/Shared/WorldMercatorProjection.cs +++ b/MapProjections/Shared/WorldMercatorProjection.cs @@ -1,5 +1,4 @@ -using System; -#if WPF +#if WPF using System.Windows; #elif AVALONIA using Avalonia; @@ -36,9 +35,7 @@ namespace MapControl.Projections public override Point RelativeScale(double latitude, double longitude) { - var lat = latitude * Math.PI / 180d; - var eSinLat = MapControl.WorldMercatorProjection.Wgs84Eccentricity * Math.Sin(lat); - var k = Math.Sqrt(1d - eSinLat * eSinLat) / Math.Cos(lat); // p.44 (7-8) + var k = MapControl.WorldMercatorProjection.ScaleFactor(latitude); return new Point(k, k); }