Post
Ask Question

The Orthodromic Distance Between Two Geo-Points

Introduction

Calculation of the great-circle (orthodromic) distance between two geo-points on the Earth surface is one of the core Geographic Information System (GIS) problems. This seemingly trivial task requires quite non-trivial algorithmic solution. Indeed, should the problem pertains to the plane geometry, then Pythagorean Theorem will provide a simple solution. But the actual GIS computations are dealing with 3d-models, namely spherical Earth representation, which requires more elaborate solution. Another level of complexity relates to more accurate ellipsoidal Earth model, which is sort of "overkill" for the majority of practical application. Spherical model results in systemic error margin within 0.3% which is acceptable for most commercial-grade apps. The second one (i.e. ellipsoidal model of the Earth ) theoretically limits error margin to the fraction of mm while dramatically increasing the computational complexity. The following solution is based on the spherical Earth model, describing 3 practical algorithms written in C#, which differ by the computational performance/accuracy.

Background

Mathematically speaking, all three algos described below result in computation of a great-circle (orthodromic) distance on Earth between 2 points, though the accuracy and performance are different. They all are based on spherical model of the Earth and provide reasonably good approximation with error margin typically not exceeding couple meters within NY City boundaries. More accurate ellipsoidal Earth model and corresponding high-accuracy Vincenty’s solution exists reducing the error margin to the fraction of mm, but also substantially increasing the computational complexity beyond the reasonable level. Therefore, 3 following algorithms based on spherical Earth model has been developed and recommended for general commercial apps, having good computational performance and reasonable accuracy.

Using the code

Below you can find three algorithmic solutions pertinent to the calculation of the great-circle (orthodromic) distance between two geo-points on the Earth surface:

  1. using System;    
  2.     
  3. namespace BusNY    
  4. {    
  5.     internal enum UnitSystem { SI = 0, US = 1 }    
  6.         
  7.     internal static class GIS    
  8.     {    
  9.         #region internal: properties (read-only)    
  10.         internal static double EarthRadiusKm { get {return _radiusEarthKM;} }    
  11.         internal static double EarthRadiusMiles { get { return _radiusEarthMiles; } }    
  12.         internal static double m2km { get { return _m2km; } }    
  13.         internal static double Deg2rad { get { return _toRad; } }    
  14.         #endregion    
  15.   
  16.         #region private: const    
  17.         private const double _radiusEarthMiles = 3959;    
  18.         private const double _radiusEarthKM = 6371;    
  19.         private const double _m2km = 1.60934;    
  20.         private const double _toRad = Math.PI / 180;    
  21.         #endregion    
  22.   
  23.         #region Method 1: Haversine algo    
  24.         /// <summary>    
  25.         /// Distance between two geographic points on surface, km/miles    
  26.         /// Haversine formula to calculate    
  27.         /// great-circle (orthodromic) distance on Earth    
  28.         /// High Accuracy, Medium speed    
  29.         /// re: http://en.wikipedia.org/wiki/Haversine_formula    
  30.         /// </summary>    
  31.         /// <param name="Lat1">double: 1st point Latitude</param>    
  32.         /// <param name="Lon1">double: 1st point Longitude</param>    
  33.         /// <param name="Lat2">double: 2nd point Latitude</param>    
  34.         /// <param name="Lon2">double: 2nd point Longitude</param>    
  35.         /// <returns>double: distance, km/miles</returns>    
  36.         internal static double DistanceHaversine(double Lat1,    
  37.                                                     double Lon1,    
  38.                                                     double Lat2,    
  39.                                                     double Lon2,    
  40.                                                     UnitSystem UnitSys ){    
  41.             try {    
  42.                 double _radLat1 = Lat1 * _toRad;    
  43.                 double _radLat2 = Lat2 * _toRad;    
  44.                 double _dLatHalf = (_radLat2 - _radLat1) / 2;    
  45.                 double _dLonHalf = Math.PI * (Lon2 - Lon1) / 360;    
  46.     
  47.                 // intermediate result    
  48.                 double _a = Math.Sin(_dLatHalf);    
  49.                 _a *= _a;    
  50.     
  51.                 // intermediate result    
  52.                 double _b = Math.Sin(_dLonHalf);    
  53.                 _b *= _b * Math.Cos(_radLat1) * Math.Cos(_radLat2);    
  54.     
  55.                 // central angle, aka arc segment angular distance    
  56.                 double _centralAngle = 2 * Math.Atan2(Math.Sqrt(_a + _b), Math.Sqrt(1 - _a - _b));    
  57.     
  58.                 // great-circle (orthodromic) distance on Earth between 2 points    
  59.                 if (UnitSys == UnitSystem.SI)  { return _radiusEarthKM * _centralAngle; }    
  60.                 else { return _radiusEarthMiles * _centralAngle; }    
  61.             }    
  62.             catch { throw; }    
  63.         }    
  64.         #endregion    
  65.   
  66.         #region Method 2: Spherical Law of Cosines    
  67.         /// <summary>    
  68.         /// Distance between two geographic points on surface, km/miles    
  69.         /// Spherical Law of Cosines formula to calculate    
  70.         /// great-circle (orthodromic) distance on Earth;    
  71.         /// High Accuracy, Medium speed    
  72.         /// re: http://en.wikipedia.org/wiki/Spherical_law_of_cosines    
  73.         /// </summary>    
  74.         /// <param name="Lat1">double: 1st point Latitude</param>    
  75.         /// <param name="Lon1">double: 1st point Longitude</param>    
  76.         /// <param name="Lat2">double: 2nd point Latitude</param>    
  77.         /// <param name="Lon2">double: 2nd point Longitude</param>    
  78.         /// <returns>double: distance, km/miles</returns>    
  79.         internal static double DistanceSLC(double Lat1,    
  80.                                         double Lon1,    
  81.                                         double Lat2,    
  82.                                         double Lon2,    
  83.                                         UnitSystem UnitSys ){    
  84.             try {    
  85.                 double _radLat1 = Lat1 * _toRad;    
  86.                 double _radLat2 = Lat2 * _toRad;    
  87.                 double _radLon1 = Lon1 * _toRad;    
  88.                 double _radLon2 = Lon2 * _toRad;    
  89.     
  90.                 // central angle, aka arc segment angular distance    
  91.                 double _centralAngle = Math.Acos(Math.Sin(_radLat1) * Math.Sin(_radLat2) +    
  92.                         Math.Cos(_radLat1) * Math.Cos(_radLat2) * Math.Cos(_radLon2 - _radLon1));    
  93.     
  94.                 // great-circle (orthodromic) distance on Earth between 2 points    
  95.                 if (UnitSys == UnitSystem.SI) { return _radiusEarthKM * _centralAngle; }    
  96.                 else { return _radiusEarthMiles * _centralAngle; }    
  97.             }    
  98.             catch { throw; }    
  99.         }    
  100.         #endregion    
  101.   
  102.         #region Method 3: Spherical Earth projection    
  103.         /// <summary>    
  104.         /// Distance between two geographic points on surface, km/miles    
  105.         /// Spherical Earth projection to a plane formula (using Pythagorean Theorem)    
  106.         /// to calculate great-circle (orthodromic) distance on Earth.    
  107.         /// central angle =    
  108.         /// Sqrt((_radLat2 - _radLat1)^2 + (Cos((_radLat1 + _radLat2)/2) * (Lon2 - Lon1))^2)    
  109.         /// Medium Accuracy, Fast,    
  110.         /// relative error less than 0.1% in search area smaller than 250 miles    
  111.         /// re: http://en.wikipedia.org/wiki/Geographical_distance    
  112.         /// </summary>    
  113.         /// <param name="Lat1">double: 1st point Latitude</param>    
  114.         /// <param name="Lon1">double: 1st point Longitude</param>    
  115.         /// <param name="Lat2">double: 2nd point Latitude</param>    
  116.         /// <param name="Lon2">double: 2nd point Longitude</param>    
  117.         /// <returns>double: distance, km/miles</returns>    
  118.         public static double DistanceSEP(double Lat1,    
  119.                                         double Lon1,    
  120.                                         double Lat2,    
  121.                                         double Lon2,    
  122.                                         UnitSystem UnitSys ){    
  123.             try    
  124.             {    
  125.                 double _radLat1 = Lat1 * _toRad;    
  126.                 double _radLat2 = Lat2 * _toRad;    
  127.                 double _dLat = (_radLat2 - _radLat1);    
  128.                 double _dLon = (Lon2 - Lon1) * _toRad;    
  129.     
  130.                 double _a = (_dLon) * Math.Cos((_radLat1 + _radLat2) / 2);    
  131.     
  132.                 // central angle, aka arc segment angular distance    
  133.                 double _centralAngle = Math.Sqrt(_a * _a + _dLat * _dLat);    
  134.     
  135.                 // great-circle (orthodromic) distance on Earth between 2 points    
  136.                 if (UnitSys == UnitSystem.SI) { return _radiusEarthKM * _centralAngle; }    
  137.                 else { return _radiusEarthMiles * _centralAngle; }    
  138.             }    
  139.             catch { throw; }    
  140.         }    
  141.         #endregion    
  142.     }    
  143. }