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:

using System;

namespace BusNY
{
    internal enum UnitSystem { SI = 0, US = 1 }

    internal static class GIS
    {
        #region internal: properties (read-only)
        internal static double EarthRadiusKm { get {return _radiusEarthKM;} }
        internal static double EarthRadiusMiles { get { return _radiusEarthMiles; } }
        internal static double m2km { get { return _m2km; } }
        internal static double Deg2rad { get { return _toRad; } }
        #endregion

        #region private: const
        private const double _radiusEarthMiles = 3959;
        private const double _radiusEarthKM = 6371;
        private const double _m2km = 1.60934;
        private const double _toRad = Math.PI / 180;
        #endregion

        #region Method 1: Haversine algo
        /// <summary>
        /// Distance between two geographic points on surface, km/miles
        /// Haversine formula to calculate
        /// great-circle (orthodromic) distance on Earth
        /// High Accuracy, Medium speed
        /// re: http://en.wikipedia.org/wiki/Haversine_formula
        /// </summary>
        /// <param name="Lat1">double: 1st point Latitude</param>
        /// <param name="Lon1">double: 1st point Longitude</param>
        /// <param name="Lat2">double: 2nd point Latitude</param>
        /// <param name="Lon2">double: 2nd point Longitude</param>
        /// <returns>double: distance, km/miles</returns>
        internal static double DistanceHaversine(double Lat1,
                                                    double Lon1,
                                                    double Lat2,
                                                    double Lon2,
                                                    UnitSystem UnitSys ){
            try {
                double _radLat1 = Lat1 * _toRad;
                double _radLat2 = Lat2 * _toRad;
                double _dLatHalf = (_radLat2 - _radLat1) / 2;
                double _dLonHalf = Math.PI * (Lon2 - Lon1) / 360;

                // intermediate result
                double _a = Math.Sin(_dLatHalf);
                _a *= _a;

                // intermediate result
                double _b = Math.Sin(_dLonHalf);
                _b *= _b * Math.Cos(_radLat1) * Math.Cos(_radLat2);

                // central angle, aka arc segment angular distance
                double _centralAngle = 2 * Math.Atan2(Math.Sqrt(_a + _b), Math.Sqrt(1 - _a - _b));

                // great-circle (orthodromic) distance on Earth between 2 points
                if (UnitSys == UnitSystem.SI)  { return _radiusEarthKM * _centralAngle; }
                else { return _radiusEarthMiles * _centralAngle; }
            }
            catch { throw; }
        }
        #endregion

        #region Method 2: Spherical Law of Cosines
        /// <summary>
        /// Distance between two geographic points on surface, km/miles
        /// Spherical Law of Cosines formula to calculate
        /// great-circle (orthodromic) distance on Earth;
        /// High Accuracy, Medium speed
        /// re: http://en.wikipedia.org/wiki/Spherical_law_of_cosines
        /// </summary>
        /// <param name="Lat1">double: 1st point Latitude</param>
        /// <param name="Lon1">double: 1st point Longitude</param>
        /// <param name="Lat2">double: 2nd point Latitude</param>
        /// <param name="Lon2">double: 2nd point Longitude</param>
        /// <returns>double: distance, km/miles</returns>
        internal static double DistanceSLC(double Lat1,
                                        double Lon1,
                                        double Lat2,
                                        double Lon2,
                                        UnitSystem UnitSys ){
            try {
                double _radLat1 = Lat1 * _toRad;
                double _radLat2 = Lat2 * _toRad;
                double _radLon1 = Lon1 * _toRad;
                double _radLon2 = Lon2 * _toRad;

                // central angle, aka arc segment angular distance
                double _centralAngle = Math.Acos(Math.Sin(_radLat1) * Math.Sin(_radLat2) +
                        Math.Cos(_radLat1) * Math.Cos(_radLat2) * Math.Cos(_radLon2 - _radLon1));

                // great-circle (orthodromic) distance on Earth between 2 points
                if (UnitSys == UnitSystem.SI) { return _radiusEarthKM * _centralAngle; }
                else { return _radiusEarthMiles * _centralAngle; }
            }
            catch { throw; }
        }
        #endregion

        #region Method 3: Spherical Earth projection
        /// <summary>
        /// Distance between two geographic points on surface, km/miles
        /// Spherical Earth projection to a plane formula (using Pythagorean Theorem)
        /// to calculate great-circle (orthodromic) distance on Earth.
        /// central angle =
        /// Sqrt((_radLat2 - _radLat1)^2 + (Cos((_radLat1 + _radLat2)/2) * (Lon2 - Lon1))^2)
        /// Medium Accuracy, Fast,
        /// relative error less than 0.1% in search area smaller than 250 miles
        /// re: http://en.wikipedia.org/wiki/Geographical_distance
        /// </summary>
        /// <param name="Lat1">double: 1st point Latitude</param>
        /// <param name="Lon1">double: 1st point Longitude</param>
        /// <param name="Lat2">double: 2nd point Latitude</param>
        /// <param name="Lon2">double: 2nd point Longitude</param>
        /// <returns>double: distance, km/miles</returns>
        public static double DistanceSEP(double Lat1,
                                        double Lon1,
                                        double Lat2,
                                        double Lon2,
                                        UnitSystem UnitSys ){
            try
            {
                double _radLat1 = Lat1 * _toRad;
                double _radLat2 = Lat2 * _toRad;
                double _dLat = (_radLat2 - _radLat1);
                double _dLon = (Lon2 - Lon1) * _toRad;

                double _a = (_dLon) * Math.Cos((_radLat1 + _radLat2) / 2);

                // central angle, aka arc segment angular distance
                double _centralAngle = Math.Sqrt(_a * _a + _dLat * _dLat);

                // great-circle (orthodromic) distance on Earth between 2 points
                if (UnitSys == UnitSystem.SI) { return _radiusEarthKM * _centralAngle; }
                else { return _radiusEarthMiles * _centralAngle; }
            }
            catch { throw; }
        }
        #endregion
    }
}