using System;
using ESRI.ArcGIS.Client.Geometry;
namespace GISProject.Extensions
{
/// <summary>
/// Extension methods for geodesic calculations.
/// </summary>
public static class Geodesic
{
private const double EarthRadius = 6378.137; //kilometers. Change to miles to return all values in miles instead
/// <summary>
/// Gets the distance between two points in Kilometers.
/// </summary>
/// <param name="start">The start point.</param>
/// <param name="end">The end point.</param>
/// <returns></returns>
public static double GetSphericalDistance(this MapPoint start, MapPoint end)
{
double lon1 = start.X / 180 * Math.PI;
double lon2 = end.X / 180 * Math.PI;
double lat1 = start.Y / 180 * Math.PI;
double lat2 = end.Y / 180 * Math.PI;
return 2 * Math.Asin(Math.Sqrt(Math.Pow((Math.Sin((lat1 - lat2) / 2)), 2) +
Math.Cos(lat1) * Math.Cos(lat2) * Math.Pow(Math.Sin((lon1 - lon2) / 2), 2))) * EarthRadius;
}
/// <summary>
/// Returns a polygon with a constant distance from the center point measured on the sphere.
/// </summary>
/// <param name="center">The center.</param>
/// <param name="distKM">Radius in kilometers.</param>
/// <returns></returns>
public static Polygon GetRadiusAsPolygon(this MapPoint center, double distKM)
{
Polyline line = GetRadius(center, distKM);
Polygon poly = new Polygon();
if (line.Paths.Count > 1)
{
PointCollection ring = line.Paths[0];
MapPoint last = ring[ring.Count - 1];
for (int i = 1; i < line.Paths.Count; i++)
{
PointCollection pnts = line.Paths[i];
ring.Add(new MapPoint(180 * Math.Sign(last.X), 90 * Math.Sign(center.Y)));
last = pnts[0];
ring.Add(new MapPoint(180 * Math.Sign(last.X), 90 * Math.Sign(center.Y)));
foreach (MapPoint p in pnts)
ring.Add(p);
last = pnts[pnts.Count - 1];
}
poly.Rings.Add(ring);
//pnts.Add(first);
}
else
{
poly.Rings.Add(line.Paths[0]);
}
if (distKM > EarthRadius * Math.PI / 2 && line.Paths.Count != 2)
{
PointCollection pnts = new PointCollection();
pnts.Add(new MapPoint(-180, -90));
pnts.Add(new MapPoint(180, -90));
pnts.Add(new MapPoint(180, 90));
pnts.Add(new MapPoint(-180, 90));
pnts.Add(new MapPoint(-180, -90));
poly.Rings.Add(pnts); //Exterior
}
return poly;
}
/// <summary>
/// Returns a polyline with a constant distance from the center point measured on the sphere.
/// </summary>
/// <param name="center">The center.</param>
/// <param name="distKM">Radius in kilometers.</param>
// <returns></returns>
public static Polyline GetRadius(this MapPoint center, double distKM)
{
Polyline line = new Polyline();
PointCollection pnts = new PointCollection();
line.Paths.Add(pnts);
for (int i = 0; i < 360; i++)
{
//double angle = i / 180.0 * Math.PI;
MapPoint p = GetPointFromHeading(center, distKM, i);
if (pnts.Count > 0)
{
MapPoint lastPoint = pnts[pnts.Count - 1];
int sign = Math.Sign(p.X);
if (Math.Abs(p.X - lastPoint.X) > 180)
{ //We crossed the date line
double lat = LatitudeAtLongitude(lastPoint, p, sign * -180);
pnts.Add(new MapPoint(sign * -180, lat));
pnts = new PointCollection();
line.Paths.Add(pnts);
pnts.Add(new MapPoint(sign * 180, lat));
}
}
pnts.Add(p);
}
pnts.Add(line.Paths[0][0]);
return line;
}
/// <summary>
/// Gets the shortest path line between two points. THe line will be following the great
/// circle described by the two points.
/// </summary>
/// <param name="start">The start point.</param>
/// <param name="end">The end point.</param>
/// <returns></returns>
public static Polyline GetGeodesicLine(this MapPoint start, MapPoint end)
{
Polyline line = new Polyline();
if (Math.Abs(end.X - start.X) <= 180) // Doesn't cross dateline
{
PointCollection pnts = GetGeodesicPoints(start, end);
line.Paths.Add(pnts);
}
else
{
double lon1 = start.X / 180 * Math.PI;
double lon2 = end.X / 180 * Math.PI;
double lat1 = start.Y / 180 * Math.PI;
double lat2 = end.Y / 180 * Math.PI;
double latA = LatitudeAtLongitude(lat1, lon1, lat2, lon2, Math.PI) / Math.PI * 180;
//double latB = LatitudeAtLongitude(lat1, lon1, lat2, lon2, -180) / Math.PI * 180;
line.Paths.Add(GetGeodesicPoints(start, new MapPoint(start.X < 0 ? -180 : 180, latA)));
line.Paths.Add(GetGeodesicPoints(new MapPoint(start.X < 0 ? 180 : -180, latA), end));
}
return line;
}
private static PointCollection GetGeodesicPoints(MapPoint start, MapPoint end)
{
double lon1 = start.X / 180 * Math.PI;
double lon2 = end.X / 180 * Math.PI;
double lat1 = start.Y / 180 * Math.PI;
double lat2 = end.Y / 180 * Math.PI;
double dX = end.X - start.X;
int points = (int)Math.Floor(Math.Abs(dX));
dX = lon2 - lon1;
PointCollection pnts = new PointCollection();
pnts.Add(start);
for (int i = 1; i < points; i++)
{
double lon = lon1 + dX / points * i;
double lat = LatitudeAtLongitude(lat1, lon1, lat2, lon2, lon);
pnts.Add(new MapPoint(lon / Math.PI * 180, lat / Math.PI * 180));
}
pnts.Add(end);
return pnts;
}
/// <summary>
/// Gets the latitude at a specific longitude for a great circle defined by p1 and p2.
/// </summary>
/// <param name="p1">The p1.</param>
/// <param name="p2">The p2.</param>
/// <param name="lon">The longitude in degrees.</param>
/// <returns></returns>
private static double LatitudeAtLongitude(MapPoint p1, MapPoint p2, double lon)
{
double lon1 = p1.X / 180 * Math.PI;
double lon2 = p2.X / 180 * Math.PI;
double lat1 = p1.Y / 180 * Math.PI;
double lat2 = p2.Y / 180 * Math.PI;
lon = lon / 180 * Math.PI;
return LatitudeAtLongitude(lat1, lon1, lat2, lon2, lon) / Math.PI * 180;
}
/// <summary>
/// Gets the latitude at a specific longitude for a great circle defined by lat1,lon1 and lat2,lon2.
/// </summary>
/// <param name="lat1">The start latitude in radians.</param>
/// <param name="lon1">The start longitude in radians.</param>
/// <param name="lat2">The end latitude in radians.</param>
/// <param name="lon2">The end longitude in radians.</param>
/// <param name="lon">The longitude in radians for where the latitude is.</param>
/// <returns></returns>
private static double LatitudeAtLongitude(double lat1, double lon1, double lat2, double lon2, double lon)
{
return Math.Atan((Math.Sin(lat1) * Math.Cos(lat2) * Math.Sin(lon - lon2)
- Math.Sin(lat2) * Math.Cos(lat1) * Math.Sin(lon - lon1)) / (Math.Cos(lat1) * Math.Cos(lat2) * Math.Sin(lon1 - lon2)));
}
/// <summary>
/// Gets the true bearing at a distance from the start point towards the new point.
/// </summary>
/// <param name="start">The start point.</param>
/// <param name="end">The point to get the bearing towards.</param>
/// <param name="distanceKM">The distance in kilometers travelled between start and end.</param>
/// <returns></returns>
public static double GetTrueBearing(MapPoint start, MapPoint end, double distanceKM)
{
double d = distanceKM / EarthRadius; //Angular distance in radians
double lon1 = start.X / 180 * Math.PI;
double lat1 = start.Y / 180 * Math.PI;
double lon2 = end.X / 180 * Math.PI;
double lat2 = end.Y / 180 * Math.PI;
double tc1;
if (Math.Sin(lon2 - lon1) < 0)
tc1 = Math.Acos((Math.Sin(lat2) - Math.Sin(lat1) * Math.Cos(d)) / (Math.Sin(d) * Math.Cos(lat1)));
else
tc1 = 2 * Math.PI - Math.Acos((Math.Sin(lat2) - Math.Sin(lat1) * Math.Cos(d)) / (Math.Sin(d) * Math.Cos(lat1)));
return tc1 / Math.PI * 180;
}
/// <summary>
/// Gets the point based on a start point, a heading and a distance.
/// </summary>
/// <param name="start">The start.</param>
/// <param name="distanceKM">The distance KM.</param>
/// <param name="heading">The heading.</param>
/// <returns></returns>
public static MapPoint GetPointFromHeading(MapPoint start, double distanceKM, double heading)
{
double brng = heading / 180 * Math.PI;
double lon1 = start.X / 180 * Math.PI;
double lat1 = start.Y / 180 * Math.PI;
double dR = distanceKM / 6378.137; //Angular distance in radians
double lat2 = Math.Asin(Math.Sin(lat1) * Math.Cos(dR) + Math.Cos(lat1) * Math.Sin(dR) * Math.Cos(brng));
double lon2 = lon1 + Math.Atan2(Math.Sin(brng) * Math.Sin(dR) * Math.Cos(lat1), Math.Cos(dR) - Math.Sin(lat1) * Math.Sin(lat2));
double lon = lon2 / Math.PI * 180;
double lat = lat2 / Math.PI * 180;
while (lon < -180) lon += 360;
while (lat < -90) lat += 180;
while (lon > 180) lon -= 360;
while (lat > 90) lat -= 180;
return new MapPoint(lon, lat);
}
}
}