diff --git a/MapControl/Shared/AzimuthalEquidistantProjection.cs b/MapControl/Shared/AzimuthalEquidistantProjection.cs index 4002e5be..279769b9 100644 --- a/MapControl/Shared/AzimuthalEquidistantProjection.cs +++ b/MapControl/Shared/AzimuthalEquidistantProjection.cs @@ -29,14 +29,15 @@ namespace MapControl public override Matrix RelativeScale(double latitude, double longitude) { var p = GetProjectedPoint(latitude, longitude); - var k = 1d; - if (p.CosC < 1d) + if (Math.Abs(p.CosC) >= 1d) { - var c = Math.Acos(p.CosC); - k = c / Math.Sin(c); // p.195 (25-2) + return new Matrix(1d, 0d, 0d, 1d, 0d, 0d); } + var c = Math.Acos(p.CosC); + var k = c / Math.Sin(c); // p.195 (25-2) + var scale = new Matrix(1d, 0d, 0d, k, 0d, 0d); // h == 1 scale.Rotate(-Math.Atan2(p.Y, p.X) * 180d / Math.PI); return scale; @@ -45,14 +46,15 @@ namespace MapControl public override Point? LocationToMap(double latitude, double longitude) { var p = GetProjectedPoint(latitude, longitude); - var k = 1d; - if (p.CosC < 1d) + if (Math.Abs(p.CosC) >= 1d) { - var c = Math.Acos(p.CosC); - k = c / Math.Sin(c); // p.195 (25-2) + return new Point(); } + var c = Math.Acos(p.CosC); + var k = c / Math.Sin(c); // p.195 (25-2) + return new Point(EarthRadius * k * p.X, EarthRadius * k * p.Y); // p.195 (22-4/5) } diff --git a/MapControl/Shared/AzimuthalProjection.cs b/MapControl/Shared/AzimuthalProjection.cs index 6f4c9412..e31d8874 100644 --- a/MapControl/Shared/AzimuthalProjection.cs +++ b/MapControl/Shared/AzimuthalProjection.cs @@ -46,14 +46,39 @@ namespace MapControl protected Location GetLocation(double x, double y, double rho, double sinC) { - var cosC = Math.Sqrt(1d - sinC * sinC); + var cos2C = 1d - sinC * sinC; + + if (cos2C < 0d) + { + return null; + } + + var cosC = Math.Sqrt(cos2C); var phi1 = Center.Latitude * 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 dLambda = Math.Atan2(x * sinC, rho * cosPhi1 * cosC - y * sinPhi1 * sinC); // (20-15), λ - λ0 + double u, v; - return new Location(phi * 180d / Math.PI, dLambda * 180d / Math.PI + Center.Longitude); + if (Center.Latitude == 90d) // (20-16) + { + u = x; + v = -y; + } + else if (Center.Latitude == -90d) // (20-17) + { + u = x; + v = y; + } + else // (20-15) + { + u = x * sinC; + v = rho * cosPhi1 * cosC - y * sinPhi1 * sinC; + } + + return new Location( + phi * 180d / Math.PI, + Math.Atan2(u, v) * 180d / Math.PI + Center.Longitude); } } } diff --git a/MapControl/Shared/GnomonicProjection.cs b/MapControl/Shared/GnomonicProjection.cs index 0f264b3e..17917f6e 100644 --- a/MapControl/Shared/GnomonicProjection.cs +++ b/MapControl/Shared/GnomonicProjection.cs @@ -41,7 +41,7 @@ namespace MapControl { var p = GetProjectedPoint(latitude, longitude); - if (p.CosC <= 0d) + if (p.CosC <= 0d) // p.167 "If cos c is zero or negative, the point is to be rejected." { return null; } diff --git a/MapControl/Shared/MapProjectionFactory.cs b/MapControl/Shared/MapProjectionFactory.cs index 645a3d6a..db4214bd 100644 --- a/MapControl/Shared/MapProjectionFactory.cs +++ b/MapControl/Shared/MapProjectionFactory.cs @@ -6,7 +6,7 @@ namespace MapControl { public virtual MapProjection GetProjection(string crsId) { - MapProjection projection = crsId switch + var projection = crsId switch { WebMercatorProjection.DefaultCrsId => new WebMercatorProjection(), WorldMercatorProjection.DefaultCrsId => new WorldMercatorProjection(), @@ -26,17 +26,27 @@ namespace MapControl return projection ?? throw new NotSupportedException($"MapProjection \"{crsId}\" is not supported."); } - public MapProjection GetProjectionFromEpsgCode(string crsId) => - crsId.StartsWith("EPSG:") && int.TryParse(crsId.Substring(5), out int epsgCode) ? GetProjection(epsgCode) : null; - - public virtual MapProjection GetProjection(int epsgCode) => epsgCode switch + public virtual MapProjection GetProjection(int epsgCode) { - var code when code >= Etrs89UtmProjection.FirstZoneEpsgCode && code <= Etrs89UtmProjection.LastZoneEpsgCode => new Etrs89UtmProjection(epsgCode % 100), - var code when code >= Nad27UtmProjection.FirstZoneEpsgCode && code <= Nad27UtmProjection.LastZoneEpsgCode => new Nad27UtmProjection(epsgCode % 100), - var code when code >= Nad83UtmProjection.FirstZoneEpsgCode && code <= Nad83UtmProjection.LastZoneEpsgCode => new Nad83UtmProjection(epsgCode % 100), - var code when code >= Wgs84UtmProjection.FirstZoneNorthEpsgCode && code <= Wgs84UtmProjection.LastZoneNorthEpsgCode => new Wgs84UtmProjection(epsgCode % 100, Hemisphere.North), - var code when code >= Wgs84UtmProjection.FirstZoneSouthEpsgCode && code <= Wgs84UtmProjection.LastZoneSouthEpsgCode => new Wgs84UtmProjection(epsgCode % 100, Hemisphere.South), - _ => null - }; + return epsgCode switch + { + var c when c is >= Etrs89UtmProjection.FirstZoneEpsgCode + and <= Etrs89UtmProjection.LastZoneEpsgCode => new Etrs89UtmProjection(c % 100), + var c when c is >= Nad27UtmProjection.FirstZoneEpsgCode + and <= Nad27UtmProjection.LastZoneEpsgCode => new Nad27UtmProjection(c % 100), + var c when c is >= Nad83UtmProjection.FirstZoneEpsgCode + and <= Nad83UtmProjection.LastZoneEpsgCode => new Nad83UtmProjection(c % 100), + var c when c is >= Wgs84UtmProjection.FirstZoneNorthEpsgCode + and <= Wgs84UtmProjection.LastZoneNorthEpsgCode => new Wgs84UtmProjection(c % 100, Hemisphere.North), + var c when c is >= Wgs84UtmProjection.FirstZoneSouthEpsgCode + and <= Wgs84UtmProjection.LastZoneSouthEpsgCode => new Wgs84UtmProjection(c % 100, Hemisphere.South), + _ => null + }; + } + + protected MapProjection GetProjectionFromEpsgCode(string crsId) + { + return crsId.StartsWith("EPSG:") && int.TryParse(crsId.Substring(5), out int epsgCode) ? GetProjection(epsgCode) : null; + } } } diff --git a/MapControl/Shared/OrthographicProjection.cs b/MapControl/Shared/OrthographicProjection.cs index df2b50d3..a7de14e6 100644 --- a/MapControl/Shared/OrthographicProjection.cs +++ b/MapControl/Shared/OrthographicProjection.cs @@ -40,7 +40,12 @@ namespace MapControl { var p = GetProjectedPoint(latitude, longitude); - return p.CosC >= 0d ? new Point(EarthRadius * p.X, EarthRadius * p.Y) : null; // p.149 (20-3/4) + if (p.CosC < 0d) // p.149 "If cos c is negative, the point is not to be plotted." + { + return null; + } + + return new Point(EarthRadius * p.X, EarthRadius * p.Y); // p.149 (20-3/4) } public override Location MapToLocation(double x, double y) @@ -48,7 +53,7 @@ namespace MapControl var rho = Math.Sqrt(x * x + y * y); var sinC = rho / EarthRadius; // p.150 (20-19) - return sinC <= 1d ? GetLocation(x, y, rho, sinC) : null; + return GetLocation(x, y, rho, sinC); } } } diff --git a/MapProjections/Shared/ProjNetMapProjectionFactory.cs b/MapProjections/Shared/ProjNetMapProjectionFactory.cs index 77ff6f31..0c130b2d 100644 --- a/MapProjections/Shared/ProjNetMapProjectionFactory.cs +++ b/MapProjections/Shared/ProjNetMapProjectionFactory.cs @@ -27,7 +27,7 @@ namespace MapControl.Projections public override MapProjection GetProjection(string crsId) { - MapProjection projection = crsId switch + return crsId switch { MapControl.WebMercatorProjection.DefaultCrsId => new WebMercatorProjection(), MapControl.WorldMercatorProjection.DefaultCrsId => new WorldMercatorProjection(), @@ -36,21 +36,33 @@ namespace MapControl.Projections MapControl.Wgs84AutoUtmProjection.DefaultCrsId => new Wgs84AutoUtmProjection(), MapControl.OrthographicProjection.DefaultCrsId => new Wgs84OrthographicProjection(), MapControl.StereographicProjection.DefaultCrsId => new Wgs84StereographicProjection(), - _ => GetProjectionFromEpsgCode(crsId) + _ => GetProjectionFromEpsgCode(crsId) ?? base.GetProjection(crsId) }; - - return projection ?? base.GetProjection(crsId); } - public override MapProjection GetProjection(int epsgCode) => epsgCode switch + public override MapProjection GetProjection(int epsgCode) { - var code when code >= Ed50UtmProjection.FirstZoneEpsgCode && code <= Ed50UtmProjection.LastZoneEpsgCode => new Ed50UtmProjection(epsgCode % 100), - var code when code >= Etrs89UtmProjection.FirstZoneEpsgCode && code <= Etrs89UtmProjection.LastZoneEpsgCode => new Etrs89UtmProjection(epsgCode % 100), - var code when code >= Nad27UtmProjection.FirstZoneEpsgCode && code <= Nad27UtmProjection.LastZoneEpsgCode => new Nad27UtmProjection(epsgCode % 100), - var code when code >= Nad83UtmProjection.FirstZoneEpsgCode && code <= Nad83UtmProjection.LastZoneEpsgCode => new Nad83UtmProjection(epsgCode % 100), - var code when code >= Wgs84UtmProjection.FirstZoneNorthEpsgCode && code <= Wgs84UtmProjection.LastZoneNorthEpsgCode => new Wgs84UtmProjection(epsgCode % 100, Hemisphere.North), - var code when code >= Wgs84UtmProjection.FirstZoneSouthEpsgCode && code <= Wgs84UtmProjection.LastZoneSouthEpsgCode => new Wgs84UtmProjection(epsgCode % 100, Hemisphere.South), - _ => CoordinateSystemWkts.TryGetValue(epsgCode, out string wkt) ? new ProjNetMapProjection(wkt) : base.GetProjection(epsgCode) - }; + if (CoordinateSystemWkts.TryGetValue(epsgCode, out string wkt)) + { + return new ProjNetMapProjection(wkt); + } + + return epsgCode switch + { + var c when c is >= Ed50UtmProjection.FirstZoneEpsgCode + and <= Ed50UtmProjection.LastZoneEpsgCode => new Ed50UtmProjection(c % 100), + var c when c is >= Etrs89UtmProjection.FirstZoneEpsgCode + and <= Etrs89UtmProjection.LastZoneEpsgCode => new Etrs89UtmProjection(c % 100), + var c when c is >= Nad27UtmProjection.FirstZoneEpsgCode + and <= Nad27UtmProjection.LastZoneEpsgCode => new Nad27UtmProjection(c % 100), + var c when c is >= Nad83UtmProjection.FirstZoneEpsgCode + and <= Nad83UtmProjection.LastZoneEpsgCode => new Nad83UtmProjection(c % 100), + var c when c is >= Wgs84UtmProjection.FirstZoneNorthEpsgCode + and <= Wgs84UtmProjection.LastZoneNorthEpsgCode => new Wgs84UtmProjection(c % 100, Hemisphere.North), + var c when c is >= Wgs84UtmProjection.FirstZoneSouthEpsgCode + and <= Wgs84UtmProjection.LastZoneSouthEpsgCode => new Wgs84UtmProjection(c % 100, Hemisphere.South), + _ => base.GetProjection(epsgCode) + }; + } } }