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);
- }
- }
-}