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;
  3. namespace BusNY
  4. {
  5. internal enum UnitSystem { SI = 0, US = 1 }
  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
  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
  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;
  47. // intermediate result
  48. double _a = Math.Sin(_dLatHalf);
  49. _a *= _a;
  51. // intermediate result
  52. double _b = Math.Sin(_dLonHalf);
  53. _b *= _b * Math.Cos(_radLat1) * Math.Cos(_radLat2);
  55. // central angle, aka arc segment angular distance
  56. double _centralAngle = 2 * Math.Atan2(Math.Sqrt(_a + _b), Math.Sqrt(1 - _a - _b));
  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
  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;
  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));
  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
  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;
  130. double _a = (_dLon) * Math.Cos((_radLat1 + _radLat2) / 2);
  132. // central angle, aka arc segment angular distance
  133. double _centralAngle = Math.Sqrt(_a * _a + _dLat * _dLat);
  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. }