Calculate distance between 2 GPS coordinates

Question

How do I calculate distance between two GPS coordinates (using latitude and longitude)?

Solution 1

Calculate the distance between two coordinates by latitude and longitude, including a Javascript implementation.

West and South locations are negative. Remember minutes and seconds are out of 60 so S31 30' is -31.50 degrees.

Don't forget to convert degrees to radians. Many languages have this function. Or its a simple calculation: radians = degrees * PI / 180.

function degreesToRadians(degrees) {
  return degrees * Math.PI / 180;
}

function distanceInKmBetweenEarthCoordinates(lat1, lon1, lat2, lon2) {
  var earthRadiusKm = 6371;

  var dLat = degreesToRadians(lat2-lat1);
  var dLon = degreesToRadians(lon2-lon1);

  lat1 = degreesToRadians(lat1);
  lat2 = degreesToRadians(lat2);

  var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
          Math.sin(dLon/2) * Math.sin(dLon/2) * Math.cos(lat1) * Math.cos(lat2); 
  var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); 
  return earthRadiusKm * c;
}

Here are some examples of usage:

distanceInKmBetweenEarthCoordinates(0,0,0,0)  // Distance between same 
                                              // points should be 0
0

distanceInKmBetweenEarthCoordinates(51.5, 0, 38.8, -77.1) // From London
                                                          // to Arlington
5918.185064088764

Solution 2

Look for haversine with Google; here is my solution:

#include <math.h>
#include "haversine.h"

#define d2r (M_PI / 180.0)

//calculate haversine distance for linear distance
double haversine_km(double lat1, double long1, double lat2, double long2)
{
    double dlong = (long2 - long1) * d2r;
    double dlat = (lat2 - lat1) * d2r;
    double a = pow(sin(dlat/2.0), 2) + cos(lat1*d2r) * cos(lat2*d2r) * pow(sin(dlong/2.0), 2);
    double c = 2 * atan2(sqrt(a), sqrt(1-a));
    double d = 6367 * c;

    return d;
}

double haversine_mi(double lat1, double long1, double lat2, double long2)
{
    double dlong = (long2 - long1) * d2r;
    double dlat = (lat2 - lat1) * d2r;
    double a = pow(sin(dlat/2.0), 2) + cos(lat1*d2r) * cos(lat2*d2r) * pow(sin(dlong/2.0), 2);
    double c = 2 * atan2(sqrt(a), sqrt(1-a));
    double d = 3956 * c; 

    return d;
}

Solution 3

C# Version of Haversine

double _eQuatorialEarthRadius = 6378.1370D;
double _d2r = (Math.PI / 180D);

private int HaversineInM(double lat1, double long1, double lat2, double long2)
{
    return (int)(1000D * HaversineInKM(lat1, long1, lat2, long2));
}

private double HaversineInKM(double lat1, double long1, double lat2, double long2)
{
    double dlong = (long2 - long1) * _d2r;
    double dlat = (lat2 - lat1) * _d2r;
    double a = Math.Pow(Math.Sin(dlat / 2D), 2D) + Math.Cos(lat1 * _d2r) * Math.Cos(lat2 * _d2r) * Math.Pow(Math.Sin(dlong / 2D), 2D);
    double c = 2D * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1D - a));
    double d = _eQuatorialEarthRadius * c;

    return d;
}

Here's a .NET Fiddle of this, so you can test it out with your own Lat/Longs.

Solution 4

Java Version of Haversine Algorithm based on Roman Makarov`s reply to this thread

public class HaversineAlgorithm {

    static final double _eQuatorialEarthRadius = 6378.1370D;
    static final double _d2r = (Math.PI / 180D);

    public static int HaversineInM(double lat1, double long1, double lat2, double long2) {
        return (int) (1000D * HaversineInKM(lat1, long1, lat2, long2));
    }

    public static double HaversineInKM(double lat1, double long1, double lat2, double long2) {
        double dlong = (long2 - long1) * _d2r;
        double dlat = (lat2 - lat1) * _d2r;
        double a = Math.pow(Math.sin(dlat / 2D), 2D) + Math.cos(lat1 * _d2r) * Math.cos(lat2 * _d2r)
                * Math.pow(Math.sin(dlong / 2D), 2D);
        double c = 2D * Math.atan2(Math.sqrt(a), Math.sqrt(1D - a));
        double d = _eQuatorialEarthRadius * c;

        return d;
    }

}

Solution 5

This is very easy to do with geography type in SQL Server 2008.

SELECT geography::Point(lat1, lon1, 4326).STDistance(geography::Point(lat2, lon2, 4326))
-- computes distance in meters using eliptical model, accurate to the mm

4326 is SRID for WGS84 elipsoidal Earth model

View more solutions

Comments

[Nadeem Shukoor]

I used haversine in my projects too.

[Hlung]

You can replace (M_PI / 180.0) with 0.0174532925199433 for better performance.

[pms]

In terms of performance: one could calculate sin(dlat/2.0) only once, store it in variable a1, and instead of pow(,2) it's MUCH better to use a1*a1. The same for the other pow(,2).

[Vinney Kelly]

In case it's not obvious, the toRad() method is a customization to the Number prototype such as: Number.prototype.toRad = function() { return this * (Math.PI / 180); }; . Or, as indicated below, you can replace (Math.PI/2) with 0.0174532925199433 (...whatever precision you deem necessary) for increased performance.

[Dylan Knowles]

If anyone, specifically those of you who don't look for end of line comments, is staring at this formula and looking for a unit of distance, the unit is km. :)

[Admin]

Yeah, or just use a post-’60s compiler.

[Pure.Krome]

I've also added a checky .NET fiddle so people can easily test this out.

[Legionar]

MySQL said Something is wrong in your syntax near '' on line 8 // declare distance float;

[Przemek]

I think your function GetDistance returns value in meters

[Veda]

Why are the cosine of lat1 and lat2 calculated and not over lon1 and lon2?

[Patrick Murphy]

@VinneyKelly Small typo but replace (Math.PI/180) not (Math.PI/2), thanks for everyones help

[Dmitry]

Does anybody have Delphi or Pascal version of the algorithm?

[Christian Kreiter]

What units is this displayed in? I printed out d, but I don't know what the number is being printed. I'm trying to figure out how to convert this to miles.

[fnx]

the .Net Framework has a build in method GeoCoordinate.GetDistanceTo. The assembly System.Device has to be referenced. MSDN Article msdn.microsoft.com/en-us/library/…

[Patrick Cornelissen]

There is no need to "optimize" (M_PI / 180.0) to a constant that no one understands without context. The compiler calculates these fixed terms for you!

[Paulo Miguel Almeida]

@Radu make sure you're using it correctly and not exchanging lat/log places when passing them to any method.

[John Machin]

This is the "spherical law of cosines" calculation which is the least accurate and most error-prone method of calculation of a great circle distance.

[John Machin]

This is the "spherical law of cosines" calculation which is the least accurate and most error-prone method of calculation of a great circle distance.

[John Machin]

This is the "spherical law of cosines" calculation which is the least accurate and most error-prone method of calculation of a great circle distance

[gub]

Is this correct? GetDirection() doesn't make use of 'dlat'.

[Nic]

@ChristianKRider Look at the first line. Think about what R normally means in math, then look up relevant, Earth-related quantities to see if the numbers match.

[Nic]

@cletus Is there a faster way to compare which of two points is farther from a third than a formula that requires sqrt and atan2 (even if the former is faster than it used to be)

[instead]

Please change lat1 and lat2 to $lat1 nad $lat2.

[chapeljuice]

For imperial units (miles) you could change earthRadiusKm to be var earthRadiusMiles = 3959;, fyi.

[Tõnu Samuel]

@hlung please test such recommendations before submitting. It does not optimize anything unless -O is disabled. Your recommendation does not work. Even more, it makes previously readable code unreadable. Which is also important.

[Hlung]

@TõnuSamuel Thank you very much for your comment. I really appreciate it. It makes sense that compiler with optimization enabled (-O) can pre-calculate operations of constants, making manual collapsing useless. I will test it out when I have time.

[Jessé Catrinck]

Another result in google maps... drive.google.com/…

[Janac Meena]

I got a reasonably close answer using this formula. I based the accuracy using this website: movable-type.co.uk/scripts/latlong.html which gave me 0.07149 km whereas your formula gave me 0.07156 which is an accuracy of about 99%

[FlorianB]

And if we wanted to give variable "a" a good name it would be halfLapsAroundGlobe (so 1 is the furthest two points can be, like Nelson (New Zealand, -41.303003, 173.256136) and Mogadouro (Portugal, 41.341175, -6.709340)

[fdermishin]

Why did you use equatorial radius instead of mean Earth radius?

[Csaba Toth]

@user13044086 Good question. It's because I derived this from Paulo Miguel Almeida's Java version. Looks like the C# version is also using that distance. Other versions here have 6371, but then you have to realize that all these algorithms may not perfectly handle the Earth's geoid shape. Feel free to modify this and use 6371. If you tell me that leads to more precise values I'll change my answer.

[fdermishin]

6371.008 is commonly used because it minimizes relative error of the formula as explained in notes on page movable-type.co.uk/scripts/latlong.html#ellipsoid

[Csaba Toth]

@user13044086 Thanks for the link, I edited my answer a while ago based on that

[Nico Haase]

Please put all information to your answer instead of linking to external ressources

[seanb]

@NicoHaase Fair call, if perhaps a tad extemporaneous - was over 12 years ago, and this was a slightly different place back then.