From 746a21854a9a97692093c07a221c477b4d921f57 Mon Sep 17 00:00:00 2001 From: ClemensF Date: Sat, 13 Feb 2021 18:50:30 +0100 Subject: [PATCH] Removed LocationEx --- MapControl/Shared/Location.cs | 50 ++++++++++++++++ MapControl/Shared/LocationCollection.cs | 76 ++++++++++++++----------- MapControl/Shared/LocationEx.cs | 59 ------------------- 3 files changed, 92 insertions(+), 93 deletions(-) delete mode 100644 MapControl/Shared/LocationEx.cs diff --git a/MapControl/Shared/Location.cs b/MapControl/Shared/Location.cs index 7ae96b4b..65febbb9 100644 --- a/MapControl/Shared/Location.cs +++ b/MapControl/Shared/Location.cs @@ -99,5 +99,55 @@ namespace MapControl return longitude; } + + /// + /// Calculates the great circle distance between this and the specified Location. + /// https://en.wikipedia.org/wiki/Great_circle + /// https://en.wikipedia.org/wiki/Great-circle_distance + /// https://en.wikipedia.org/wiki/Great-circle_navigation + /// + public double GetDistance( + Location location, double earthRadius = MapProjection.Wgs84EquatorialRadius) + { + var lat1 = latitude * Math.PI / 180d; + var lon1 = longitude * Math.PI / 180d; + var lat2 = location.latitude * Math.PI / 180d; + var lon2 = location.longitude * Math.PI / 180d; + var sinLat1 = Math.Sin(lat1); + var cosLat1 = Math.Cos(lat1); + var sinLat2 = Math.Sin(lat2); + var cosLat2 = Math.Cos(lat2); + var sinLon12 = Math.Sin(lon2 - lon1); + var cosLon12 = Math.Cos(lon2 - lon1); + var a = cosLat1 * sinLat2 - sinLat1 * cosLat2 * cosLon12; + var b = cosLat2 * sinLon12; + var s12 = Math.Atan2(Math.Sqrt(a * a + b * b), sinLat1 * sinLat2 + cosLat1 * cosLat2 * cosLon12); + + return earthRadius * s12; + } + + /// + /// Calculates the Location on a great circle at the specified azimuth angle and distance from this Location. + /// https://en.wikipedia.org/wiki/Great_circle + /// https://en.wikipedia.org/wiki/Great-circle_navigation + /// + public Location GetLocation( + double azimuth, double distance, double earthRadius = MapProjection.Wgs84EquatorialRadius) + { + var s12 = distance / earthRadius; + var az1 = azimuth * Math.PI / 180d; + var lat1 = latitude * Math.PI / 180d; + var lon1 = longitude * Math.PI / 180d; + var sinS12 = Math.Sin(s12); + var cosS12 = Math.Cos(s12); + var sinAz1 = Math.Sin(az1); + var cosAz1 = Math.Cos(az1); + var sinLat1 = Math.Sin(lat1); + var cosLat1 = Math.Cos(lat1); + var lat2 = Math.Asin(sinLat1 * cosS12 + cosLat1 * sinS12 * cosAz1); + var lon2 = lon1 + Math.Atan2(sinS12 * sinAz1, (cosLat1 * cosS12 - sinLat1 * sinS12 * cosAz1)); + + return new Location(lat2 * 180d / Math.PI, lon2 * 180d / Math.PI); + } } } diff --git a/MapControl/Shared/LocationCollection.cs b/MapControl/Shared/LocationCollection.cs index 829a81b7..020e198b 100644 --- a/MapControl/Shared/LocationCollection.cs +++ b/MapControl/Shared/LocationCollection.cs @@ -58,9 +58,12 @@ namespace MapControl } /// - /// see https://en.wikipedia.org/wiki/Great-circle_navigation + /// Calculates a series of Locations on a great circle, or orthodrome, that connects the two specified Locations, + /// with an optional angular resolution specified in degrees. + /// + /// See https://en.wikipedia.org/wiki/Great-circle_navigation /// - public static LocationCollection CalculateGreatCircleLocations(Location location1, Location location2, double resolution = 1d) + public static LocationCollection OrthodromeLocations(Location location1, Location location2, double resolution = 1d) { if (resolution <= 0d) { @@ -72,6 +75,7 @@ namespace MapControl var lon1 = location1.Longitude * Math.PI / 180d; var lat2 = location2.Latitude * Math.PI / 180d; var lon2 = location2.Longitude * Math.PI / 180d; + var cosLat1 = Math.Cos(lat1); var sinLat1 = Math.Sin(lat1); var cosLat2 = Math.Cos(lat2); @@ -83,14 +87,12 @@ namespace MapControl var b = cosLat2 * sinLon12; var s12 = Math.Atan2(Math.Sqrt(a * a + b * b), sinLat1 * sinLat2 + cosLat1 * cosLat2 * cosLon12); + var n = (int)Math.Ceiling(s12 / resolution * 180d / Math.PI); // s12 in radians + var locations = new LocationCollection(new Location(location1.Latitude, location1.Longitude)); - resolution *= Math.PI / 180d; - - if (s12 > resolution) // s12 and resolution in radians + if (n > 1) { - var n = (int)Math.Round(s12 / resolution); - var az1 = Math.Atan2(sinLon12, cosLat1 * sinLat2 / cosLat2 - sinLat1 * cosLon12); var cosAz1 = Math.Cos(az1); var sinAz1 = Math.Sin(az1); @@ -120,9 +122,12 @@ namespace MapControl } /// - /// see https://en.wikipedia.org/wiki/Rhumb_line + /// Calculates a series of Locations on a rhumb line, or loxodrome, that connects the two specified Locations, + /// with an optional angular resolution specified in degrees. + /// + /// See https://en.wikipedia.org/wiki/Rhumb_line /// - public static LocationCollection CalculateRhumbLineLocations(Location location1, Location location2, double resolution = 1d) + public static LocationCollection LoxodromeLocations(Location location1, Location location2, double resolution = 1d) { if (resolution <= 0d) { @@ -131,7 +136,10 @@ namespace MapControl } var lat1 = location1.Latitude; + var lon1 = location1.Longitude; var lat2 = location2.Latitude; + var lon2 = location2.Longitude; + var y1 = WebMercatorProjection.LatitudeToY(lat1); var y2 = WebMercatorProjection.LatitudeToY(lat2); @@ -147,50 +155,50 @@ namespace MapControl nameof(location2), "The location2 argument must have an absolute latitude value of less than 90."); } - var lon1 = location1.Longitude; - var lon2 = location2.Longitude; var dlat = lat2 - lat1; var dlon = lon2 - lon1; var dy = y2 - y1; - double s12; + // beta = atan(dlon,dy) // sec(beta) = 1 / cos(atan(dlon,dy)) = sqrt(1 + (dlon/dy)^2) + var sec = Math.Sqrt(1d + dlon * dlon / (dy * dy)); - if (sec > 1000d) // beta near +/-90° - { - var lat = (lat1 + lat2) * Math.PI / 360d; + const double secLimit = 1000d; // beta approximately +/-90° - s12 = Math.Abs(dlon * Math.Cos(lat)); + double s12; + + if (sec > secLimit) + { + var lat = (lat1 + lat2) * Math.PI / 360d; // mean latitude + + s12 = Math.Abs(dlon * Math.Cos(lat)); // distance in degrees along parallel of latitude } else { - s12 = Math.Abs(dlat * sec); + s12 = Math.Abs(dlat * sec); // distance in degrees along loxodrome } + var n = (int)Math.Ceiling(s12 / resolution); + var locations = new LocationCollection(new Location(lat1, lon1)); - if (s12 > resolution) // s12 and resolution in degress + if (sec > secLimit) { - var n = (int)Math.Round(s12 / resolution); - - if (sec > 1000d) + for (var i = 1; i < n; i++) { - for (var i = 1; i < n; i++) - { - var lon = lon1 + i * dlon / n; - var lat = WebMercatorProjection.YToLatitude(y1 + i * dy / n); - locations.Add(lat, lon); - } + var lon = lon1 + i * dlon / n; + var lat = WebMercatorProjection.YToLatitude(y1 + i * dy / n); + locations.Add(lat, lon); } - else + } + else + { + for (var i = 1; i < n; i++) { - for (var i = 1; i < n; i++) - { - var lat = lat1 + i * dlat / n; - var lon = lon1 + dlon * (WebMercatorProjection.LatitudeToY(lat) - y1) / dy; - locations.Add(lat, lon); - } + var lat = lat1 + i * dlat / n; + var lon = lon1 + dlon * (WebMercatorProjection.LatitudeToY(lat) - y1) / dy; + locations.Add(lat, lon); } } diff --git a/MapControl/Shared/LocationEx.cs b/MapControl/Shared/LocationEx.cs deleted file mode 100644 index 57e35352..00000000 --- a/MapControl/Shared/LocationEx.cs +++ /dev/null @@ -1,59 +0,0 @@ -// XAML Map Control - https://github.com/ClemensFischer/XAML-Map-Control -// © 2021 Clemens Fischer -// Licensed under the Microsoft Public License (Ms-PL) - -using System; - -namespace MapControl -{ - /// - /// Provides helper methods for geodetic calculations on a sphere. - /// - public static class LocationEx - { - /// - /// see https://en.wikipedia.org/wiki/Great-circle_navigation - /// - public static double GreatCircleDistance( - this Location location1, Location location2, double earthRadius = MapProjection.Wgs84EquatorialRadius) - { - var lat1 = location1.Latitude * Math.PI / 180d; - var lon1 = location1.Longitude * Math.PI / 180d; - var lat2 = location2.Latitude * Math.PI / 180d; - var lon2 = location2.Longitude * Math.PI / 180d; - var sinLat1 = Math.Sin(lat1); - var cosLat1 = Math.Cos(lat1); - var sinLat2 = Math.Sin(lat2); - var cosLat2 = Math.Cos(lat2); - var sinLon12 = Math.Sin(lon2 - lon1); - var cosLon12 = Math.Cos(lon2 - lon1); - var a = cosLat1 * sinLat2 - sinLat1 * cosLat2 * cosLon12; - var b = cosLat2 * sinLon12; - var s12 = Math.Atan2(Math.Sqrt(a * a + b * b), sinLat1 * sinLat2 + cosLat1 * cosLat2 * cosLon12); - - return earthRadius * s12; - } - - /// - /// see https://en.wikipedia.org/wiki/Great-circle_navigation - /// - public static Location GreatCircleLocation( - this Location location, double azimuth, double distance, double earthRadius = MapProjection.Wgs84EquatorialRadius) - { - var s12 = distance / earthRadius; - var az1 = azimuth * Math.PI / 180d; - var lat1 = location.Latitude * Math.PI / 180d; - var lon1 = location.Longitude * Math.PI / 180d; - var sinS12 = Math.Sin(s12); - var cosS12 = Math.Cos(s12); - var sinAz1 = Math.Sin(az1); - var cosAz1 = Math.Cos(az1); - var sinLat1 = Math.Sin(lat1); - var cosLat1 = Math.Cos(lat1); - var lat2 = Math.Asin(sinLat1 * cosS12 + cosLat1 * sinS12 * cosAz1); - var lon2 = lon1 + Math.Atan2(sinS12 * sinAz1, (cosLat1 * cosS12 - sinLat1 * sinS12 * cosAz1)); - - return new Location(lat2 * 180d / Math.PI, lon2 * 180d / Math.PI); - } - } -}