Calculating the Distance Between two GPS Coordinates

C#

Posted by Alex Peta on March 02, 2013 Copyright© from Bing images : Aeonium leaf detail (© Tim Gainey/Alamy)

Although this article is still related to the Windows Phone App series, I think it is stand alone subject.

Ok, so I had to calculate the distance between 2 GPS positions.

My first intention was to go with the basic formula that you learn in school. But it isn’t right. This is ok for a 2D space, not for GPS positions where we have to take in consideration the curvature of the Earth as well.

After searching the problem a bit,I struck gold The Harvesine Formula and its implementation : http://www.movable-type.co.uk/scripts/latlong.html – this is the best reference with explanations and sample code written in JavaScript.

TheHarvesine Formula

Haversine
formula:
a = sin²(?f/2) + cos(f1) * cos(f2) * sin²(??/2)
c = 2 * atan2(va, v(1-a))
d = R * c
where f is latitude, ? is longitude, R is earth’s radius (mean radius = 6,371km)
  note that angles need to be in radians to pass to trig functions!

All I had to do migrate this in C# and test it out.

I’m going to this the fancy way around because I have been brainwashed with OOP :)

A simple GPSPoint strunct that will hold the coordinates

namespace AlexPetaBlog.GPSDistanceCalculator
{
    public struct GPSPoint
    {
        #region Private Members
        private double _latitude;
        private double _longitude;
        #endregion Private Members

        #region Public Properties
        public double Latitude
        {
            get { return _latitude; }
            set { _latitude = value; }
        }
        public double Longitude
        {
            get { return _longitude; }
            set { _longitude = value; }
        }
        #endregion Public Properties

        #region Constructors
        public GPSPoint(double latitude, double longitude)
        {
            _latitude = latitude;
            _longitude = longitude;
        }
        #endregion Constructors

        #region NullObject Pattern
        private static GPSPoint _empty = new GPSPoint(0, 0);
        public static GPSPoint NullPoint = _empty;
        #endregion NullObject Pattern
    }
}

The GPSDistanceHelper class

namespace AlexPetaBlog.GPSDistanceCalculator
{
    /// 
    /// This uses the ‘haversine’ formula to calculate the great-circle distance between two points – that is, the shortest distance over the earth’s surface – giving an ‘as-the-crow-flies’ distance between the points (ignoring any hills, of course!).
    /// 
    public class GPSDistanceHelper
    {

        #region Private Members
        //represents the radius of the Earth in KM
        private readonly int R = 6371;
        private GPSPoint _start;
        private GPSPoint _end;
        private double _distance;
        #endregion Private Members

        #region Public Properties
        public GPSPoint Start
        {
            get { return _start; }
            set { _start = value; }
        }
        public GPSPoint End
        {
            get { return _end; }
            set { _end = value; }
        }
        public double Distance
        {
            get
            {
                return _distance;   
            }
            set
            {
                _distance = value;
            }
        }
        #endregion Public Properties

        #region Constructors
        public GPSDistanceHelper(): this(GPSPoint.NullPoint,GPSPoint.NullPoint)
        {
        }
        public GPSDistanceHelper(GPSPoint start,GPSPoint end)
        {
            Start = start;
            End = end;

            CalculateDistance();
        }
        #endregion Constructors



        #region Statics
        public static  double DegreesToRadians(double degrees)
        {
            return degrees * (Math.PI / 180);
        }
        #endregion Statics

        #region Private Methods
        private void CalculateDistance()
        {
            double dLat = GPSDistanceHelper.DegreesToRadians(End.Latitude - Start.Latitude);
            double dLon = GPSDistanceHelper.DegreesToRadians(End.Longitude - Start.Longitude);

            double lat1 = GPSDistanceHelper.DegreesToRadians(Start.Latitude);
            double lat2 = GPSDistanceHelper.DegreesToRadians(End.Latitude);

            double a = Math.Sin(dLat / 2) * Math.Sin(dLat / 2) + Math.Sin(dLon / 2) * Math.Sin(dLon / 2) * Math.Cos(lat1) * Math.Cos(lat2);

            double c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));

            this.Distance = R * c;
        }
        #endregion Private Methods
    }
}

And of course now we need to test this out with some unit tests.

namespace AlexPetaBlog.GPSDistanceCalculator.Tests
{
    [TestClass]
    public class GPSCalculatorTests
    {
        #region Private Members
        private GPSDistanceHelper _distanceHelper;

        private GPSPoint _bucharest;
        private GPSPoint _centralPark;

        #endregion Private Members

        [TestInitialize]
        public void Setup()
        {
            _distanceHelper = new GPSDistanceHelper();
            _bucharest = new GPSPoint(44.4269D, 026.1025D);
            _centralPark = new GPSPoint(40.7939D, -073.9571D);
        }

        [TestCleanup]
        public void TearDown()
        {
            _distanceHelper = null;
            _bucharest = GPSPoint.NullPoint;
            _centralPark = GPSPoint.NullPoint;
        }
        

        [TestMethod]
        public void GPSDistanceHelper_ToRadians_ReturnsCorrectResult()
        {
            //Arrange
            double expected = 0.174532925D;
            double delta = 0.000000001D;

            //Act
            double result = GPSDistanceHelper.DegreesToRadians(10);
            
            //Assert
            Assert.AreEqual(expected, result, delta);
        }

        [TestMethod]
        public void GPSDistanceHelper_Calculate_ReturnsCorrectResult()
        {
            //Arrange
            double expected = 7641D;
            double result = 0D;
            double delta = 0.9D;

            GPSDistanceHelper helper = new GPSDistanceHelper(start : _bucharest, end :_centralPark);

            //Act
            result = helper.Distance;

            //Assert
            Assert.AreEqual(expected, result,delta);
        }
    }
}

And the results :

Conclusion

You can see the full source code in action on Github.