# Calculate distance between 2 GPS coordinates

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

• This algorithm is known as the Great Circle distance. Dec 13, 2008 at 22:14
• @GregHewgill, the first sentence of that article says "This article is about shortest-distance on a sphere." Ie clearly not applicable to GPS coordinates. Dec 11, 2021 at 23:56

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);

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
``````
• 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. Jul 23, 2013 at 6:05
• 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. :) Sep 27, 2013 at 18:01
• @VinneyKelly Small typo but replace (Math.PI/180) not (Math.PI/2), thanks for everyones help Jul 21, 2015 at 16:44
• @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
Feb 12, 2017 at 13:01
• For imperial units (miles) you could change `earthRadiusKm` to be `var earthRadiusMiles = 3959;`, fyi. Aug 9, 2017 at 21:53

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;
}
``````
• You can replace (M_PI / 180.0) with 0.0174532925199433 for better performance. Aug 1, 2011 at 9:19
• 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).
– pms
Oct 27, 2011 at 23:34
• Yeah, or just use a post-’60s compiler.
– user1804599
Jan 28, 2014 at 10:45
• 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! Sep 19, 2016 at 5:35
• @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. Jan 2, 2018 at 15:47

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.

• I've also added a checky .NET fiddle so people can easily test this out. Aug 15, 2014 at 5:42
• 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/…
– fnx
Feb 14, 2016 at 18:12

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;
}

}
``````
• @Radu make sure you're using it correctly and not exchanging lat/log places when passing them to any method. Dec 26, 2016 at 17:48
• 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% Sep 18, 2019 at 20:40

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

Here's a Haversine function in Python that I use:

``````from math import pi,sqrt,sin,cos,atan2

def haversine(pos1, pos2):
lat1 = float(pos1['lat'])
long1 = float(pos1['long'])
lat2 = float(pos2['lat'])
long2 = float(pos2['long'])

degree_to_rad = float(pi / 180.0)

d_lat = (lat2 - lat1) * degree_to_rad
d_long = (long2 - long1) * degree_to_rad

a = pow(sin(d_lat / 2), 2) + cos(lat1 * degree_to_rad) * cos(lat2 * degree_to_rad) * pow(sin(d_long / 2), 2)
c = 2 * atan2(sqrt(a), sqrt(1 - a))
km = 6367 * c
mi = 3956 * c

return {"km":km, "miles":mi}
``````
• The math module contains a function named radians which converts from degrees to radians. `from math import radians` Jun 7 at 15:57

I needed to calculate a lot of distances between the points for my project, so I went ahead and tried to optimize the code, I have found here. On average in different browsers my new implementation runs 2 times faster than the most upvoted answer.

``````function distance(lat1, lon1, lat2, lon2) {
var p = 0.017453292519943295;    // Math.PI / 180
var c = Math.cos;
var a = 0.5 - c((lat2 - lat1) * p)/2 +
c(lat1 * p) * c(lat2 * p) *
(1 - c((lon2 - lon1) * p))/2;

return 12742 * Math.asin(Math.sqrt(a)); // 2 * R; R = 6371 km
}
``````

You can play with my jsPerf and see the results here.

Recently I needed to do the same in python, so here is a python implementation:

``````from math import cos, asin, sqrt
def distance(lat1, lon1, lat2, lon2):
p = 0.017453292519943295
a = 0.5 - cos((lat2 - lat1) * p)/2 + cos(lat1 * p) * cos(lat2 * p) * (1 - cos((lon2 - lon1) * p)) / 2
return 12742 * asin(sqrt(a))
``````

And for the sake of completeness: Haversine on wiki.

It depends on how accurate you need it to be. If you need pinpoint accuracy, it is best to look at an algorithm which uses an ellipsoid, rather than a sphere, such as Vincenty's algorithm, which is accurate to the mm.

• Please put all information to your answer instead of linking to external ressources Mar 15, 2021 at 14:12
• @NicoHaase Fair call, if perhaps a tad extemporaneous - was over 12 years ago, and this was a slightly different place back then. Mar 16, 2021 at 2:10

Here it is in C# (lat and long in radians):

``````double CalculateGreatCircleDistance(double lat1, double long1, double lat2, double long2, double radius)
{
return radius * Math.Acos(
Math.Sin(lat1) * Math.Sin(lat2)
+ Math.Cos(lat1) * Math.Cos(lat2) * Math.Cos(long2 - long1));
}
``````

If your lat and long are in degrees then divide by 180/PI to convert to radians.

• 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. Jan 11, 2017 at 2:27

PHP version:

(Remove all `deg2rad()` if your coordinates are already in radians.)

``````\$R = 6371; // km

\$a = sin(\$dLat/2) * sin(\$dLat/2) +
sin(\$dLon/2) * sin(\$dLon/2) * cos(\$lat1) * cos(\$lat2);

\$c = 2 * atan2(sqrt(\$a), sqrt(1-\$a));
\$d = \$R * \$c;
``````
• Please change lat1 and lat2 to \$lat1 nad \$lat2. Jul 31, 2017 at 5:56

A T-SQL function, that I use to select records by distance for a center

``````Create Function  [dbo].[DistanceInMiles]
(  @fromLatitude float ,
@fromLongitude float ,
@toLatitude float,
@toLongitude float
)
returns float
AS
BEGIN
declare @distance float

)as float)
return  round(@distance,1)
END
``````
• 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. Jan 11, 2017 at 2:30

I. Regarding "Breadcrumbs" method

1. Earth radius is different on different Lat. This must be taken into consideration in Haversine algorithm.
2. Consider Bearing change, which turns straight lines to arches (which are longer)
3. Taking Speed change into account will turn arches to spirals (which are longer or shorter than arches)
4. Altitude change will turn flat spirals to 3D spirals (which are longer again). This is very important for hilly areas.

Below see the function in C which takes #1 and #2 into account:

``````double   calcDistanceByHaversine(double rLat1, double rLon1, double rHeading1,
double rLat2, double rLon2, double rHeading2){
double rDLatRad = 0.0;
double rDLonRad = 0.0;
double rLat1Rad = 0.0;
double rLat2Rad = 0.0;
double a = 0.0;
double c = 0.0;
double rResult = 0.0;
double rEarthRadius = 0.0;
double rDHeading = 0.0;

if ((rLat1 < -90.0) || (rLat1 > 90.0) || (rLat2 < -90.0) || (rLat2 > 90.0)
|| (rLon1 < -180.0) || (rLon1 > 180.0) || (rLon2 < -180.0)
|| (rLon2 > 180.0)) {
return -1;
};

rDLatRad = (rLat2 - rLat1) * DEGREE_TO_RADIANS;
rDLonRad = (rLon2 - rLon1) * DEGREE_TO_RADIANS;

a = sin(rDLatRad / 2) * sin(rDLatRad / 2) + sin(rDLonRad / 2) * sin(

if (a == 0.0) {
return 0.0;
}

c = 2 * atan2(sqrt(a), sqrt(1 - a));
rEarthRadius = 6378.1370 - (21.3847 * 90.0 / ((fabs(rLat1) + fabs(rLat2))
/ 2.0));
rResult = rEarthRadius * c;

// Chord to Arc Correction based on Heading changes. Important for routes with many turns and U-turns

if ((rHeading1 >= 0.0) && (rHeading1 < 360.0) && (rHeading2 >= 0.0)
&& (rHeading2 < 360.0)) {
if (rDHeading > 180.0) {
}
if (rDHeading > 5.0) {
} else {
}
}
return rResult;
}
``````

II. There is an easier way which gives pretty good results.

By Average Speed.

Trip_distance = Trip_average_speed * Trip_time

Since GPS Speed is detected by Doppler effect and is not directly related to [Lon,Lat] it can be at least considered as secondary (backup or correction) if not as main distance calculation method.

If you need something more accurate then have a look at this.

Vincenty's formulae are two related iterative methods used in geodesy to calculate the distance between two points on the surface of a spheroid, developed by Thaddeus Vincenty (1975a) They are based on the assumption that the figure of the Earth is an oblate spheroid, and hence are more accurate than methods such as great-circle distance which assume a spherical Earth.

The first (direct) method computes the location of a point which is a given distance and azimuth (direction) from another point. The second (inverse) method computes the geographical distance and azimuth between two given points. They have been widely used in geodesy because they are accurate to within 0.5 mm (0.020″) on the Earth ellipsoid.

If you're using .NET don't reivent the wheel. See System.Device.Location. Credit to fnx in the comments in another answer.

``````using System.Device.Location;

double lat1 = 45.421527862548828D;
double long1 = -75.697189331054688D;
double lat2 = 53.64135D;
double long2 = -113.59273D;

GeoCoordinate geo1 = new GeoCoordinate(lat1, long1);
GeoCoordinate geo2 = new GeoCoordinate(lat2, long2);

double distance = geo1.GetDistanceTo(geo2);
``````

This is version from "Henry Vilinskiy" adapted for MySQL and Kilometers:

``````CREATE FUNCTION `CalculateDistanceInKm`(
fromLatitude float,
fromLongitude float,
toLatitude float,
toLongitude float
) RETURNS float
BEGIN
declare distance float;

select
6367 * ACOS(
round(
,15)
)
into distance;

return  round(distance,3);
END;
``````
• `MySQL` said `Something is wrong in your syntax near '' on line 8` `// declare distance float;` Dec 2, 2014 at 13:37
• 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 Jan 11, 2017 at 2:30

here is the Swift implementation from the answer

``````func degreesToRadians(degrees: Double) -> Double {
return degrees * Double.pi / 180
}

func distanceInKmBetweenEarthCoordinates(lat1: Double, lon1: Double, lat2: Double, lon2: Double) -> Double {

let earthRadiusKm: Double = 6371

let dLat = degreesToRadians(degrees: lat2 - lat1)
let dLon = degreesToRadians(degrees: lon2 - lon1)

let lat1 = degreesToRadians(degrees: lat1)
let lat2 = degreesToRadians(degrees: lat2)

let a = sin(dLat/2) * sin(dLat/2) +
sin(dLon/2) * sin(dLon/2) * cos(lat1) * cos(lat2)
let c = 2 * atan2(sqrt(a), sqrt(1 - a))
return earthRadiusKm * c
}
``````

This Lua code is adapted from stuff found on Wikipedia and in Robert Lipe's GPSbabel tool:

``````local EARTH_RAD = 6378137.0
-- earth's radius in meters (official geoid datum, not 20,000km / pi)

-- earth's radius in miles

local multipliers = {
radians = 1, miles = radmiles, mi = radmiles, feet = radmiles * 5280,
meters = EARTH_RAD, m = EARTH_RAD, km = EARTH_RAD / 1000,
degrees = 360 / (2 * math.pi), min = 60 * 360 / (2 * math.pi)
}

function gcdist(pt1, pt2, units) -- return distance in radians or given units
--- this formula works best for points close together or antipodal
--- rounding error strikes when distance is one-quarter Earth's circumference
--- (ref: wikipedia Great-circle distance)
if not pt1.radians then pt1 = rad(pt1) end
if not pt2.radians then pt2 = rad(pt2) end
local sdlat = sin((pt1.lat - pt2.lat) / 2.0);
local sdlon = sin((pt1.lon - pt2.lon) / 2.0);
local res = sqrt(sdlat * sdlat + cos(pt1.lat) * cos(pt2.lat) * sdlon * sdlon);
res = res > 1 and 1 or res < -1 and -1 or res
res = 2 * asin(res);
if units then return res * assert(multipliers[units])
else return res
end
end
``````
``````    private double deg2rad(double deg)
{
return (deg * Math.PI / 180.0);
}

{
return (rad / Math.PI * 180.0);
}

private double GetDistance(double lat1, double lon1, double lat2, double lon2)
{
//code for Distance in Kilo Meter
double theta = lon1 - lon2;
dist = Math.Abs(Math.Round(rad2deg(Math.Acos(dist)) * 60 * 1.1515 * 1.609344 * 1000, 0));
return (dist);
}

private double GetDirection(double lat1, double lon1, double lat2, double lon2)
{
//code for Direction in Degrees
double y = Math.Sin(dlon) * Math.Cos(lat2);
double direct = Math.Round(rad2deg(Math.Atan2(y, x)), 0);
if (direct < 0)
direct = direct + 360;
return (direct);
}

private double GetSpeed(double lat1, double lon1, double lat2, double lon2, DateTime CurTime, DateTime PrevTime)
{
//code for speed in Kilo Meter/Hour
TimeSpan TimeDifference = CurTime.Subtract(PrevTime);
double TimeDifferenceInSeconds = Math.Round(TimeDifference.TotalSeconds, 0);
double theta = lon1 - lon2;
dist = rad2deg(Math.Acos(dist)) * 60 * 1.1515 * 1.609344;
double Speed = Math.Abs(Math.Round((dist / Math.Abs(TimeDifferenceInSeconds)) * 60 * 60, 0));
return (Speed);
}

private double GetDuration(DateTime CurTime, DateTime PrevTime)
{
//code for speed in Kilo Meter/Hour
TimeSpan TimeDifference = CurTime.Subtract(PrevTime);
double TimeDifferenceInSeconds = Math.Abs(Math.Round(TimeDifference.TotalSeconds, 0));
return (TimeDifferenceInSeconds);
}
``````
• I think your function GetDistance returns value in meters Dec 21, 2014 at 14:11
• Is this correct? GetDirection() doesn't make use of 'dlat'.
– gub
Jan 28, 2017 at 12:52

i took the top answer and used it in a Scala program

``````import java.lang.Math.{atan2, cos, sin, sqrt}

def latLonDistance(lat1: Double, lon1: Double)(lat2: Double, lon2: Double): Double = {
val earthRadiusKm = 6371
val dLat = (lat2 - lat1).toRadians
val dLon = (lon2 - lon1).toRadians

val a = sin(dLat / 2) * sin(dLat / 2) + sin(dLon / 2) * sin(dLon / 2) * cos(latRad1) * cos(latRad2)
val c = 2 * atan2(sqrt(a), sqrt(1 - a))
}
``````

i curried the function in order to be able to easily produce functions that have one of the two locations fixed and require only a pair of lat/lon to produce distance.

Here's a Kotlin variation:

``````import kotlin.math.*

class HaversineAlgorithm {

companion object {
private const val MEAN_EARTH_RADIUS = 6371.008
private const val D2R = Math.PI / 180.0
}

private fun haversineInKm(lat1: Double, lon1: Double, lat2: Double, lon2: Double): Double {
val lonDiff = (lon2 - lon1) * D2R
val latDiff = (lat2 - lat1) * D2R
val latSin = sin(latDiff / 2.0)
val lonSin = sin(lonDiff / 2.0)
val a = latSin * latSin + (cos(lat1 * D2R) * cos(lat2 * D2R) * lonSin * lonSin)
val c = 2.0 * atan2(sqrt(a), sqrt(1.0 - a))
return MEAN_EARTH_RADIUS * c
}
}
``````
• Why did you use equatorial radius instead of mean Earth radius? May 25, 2020 at 0:47
• @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. May 25, 2020 at 6:34
• 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 May 25, 2020 at 7:20
• @user13044086 Thanks for the link, I edited my answer a while ago based on that Jun 8, 2020 at 0:45

I guess you want it along the curvature of the earth. Your two points and the center of the earth are on a plane. The center of the earth is the center of a circle on that plane and the two points are (roughly) on the perimeter of that circle. From that you can calculate the distance by finding out what the angle from one point to the other is.

If the points are not the same heights, or if you need to take into account that the earth is not a perfect sphere it gets a little more difficult.

you can find a implementation of this (with some good explanation) in F# on fssnip

here are the important parts:

``````
let GreatCircleDistance<[<Measure>] 'u> (R : float<'u>) (p1 : Location) (p2 : Location) =
let degToRad (x : float<deg>) = System.Math.PI * x / 180.0<deg/rad>

let sq x = x * x
// take the sin of the half and square the result
let sinSqHf (a : float<rad>) = (System.Math.Sin >> sq) (a / 2.0<rad>)
let cos (a : float<deg>) = System.Math.Cos (degToRad a / 1.0<rad>)

let dLat = (p2.Latitude - p1.Latitude) |> degToRad
let dLon = (p2.Longitude - p1.Longitude) |> degToRad

let a = sinSqHf dLat + cos p1.Latitude * cos p2.Latitude * sinSqHf dLon
let c = 2.0 * System.Math.Atan2(System.Math.Sqrt(a), System.Math.Sqrt(1.0-a))

R * c
``````

I needed to implement this in PowerShell, hope it can help someone else. Some notes about this method

1. Don't split any of the lines or the calculation will be wrong
2. To calculate in KM remove the * 1000 in the calculation of \$distance
3. Change \$earthsRadius = 3963.19059 and remove * 1000 in the calculation of \$distance the to calulate the distance in miles
4. I'm using Haversine, as other posts have pointed out Vincenty's formulae is much more accurate

``````Function MetresDistanceBetweenTwoGPSCoordinates(\$latitude1, \$longitude1, \$latitude2, \$longitude2)
{
\$Rad = ([math]::PI / 180);

\$earthsRadius = 6378.1370 # Earth's Radius in KM
\$dLat = (\$latitude2 - \$latitude1) * \$Rad
\$dLon = (\$longitude2 - \$longitude1) * \$Rad
\$latitude1 = \$latitude1 * \$Rad
\$latitude2 = \$latitude2 * \$Rad

\$a = [math]::Sin(\$dLat / 2) * [math]::Sin(\$dLat / 2) + [math]::Sin(\$dLon / 2) * [math]::Sin(\$dLon / 2) * [math]::Cos(\$latitude1) * [math]::Cos(\$latitude2)
\$c = 2 * [math]::ATan2([math]::Sqrt(\$a), [math]::Sqrt(1-\$a))

\$distance = [math]::Round(\$earthsRadius * \$c * 1000, 0) #Multiple by 1000 to get metres

Return \$distance
}
``````

Scala version

``````  def deg2rad(deg: Double) = deg * Math.PI / 180.0

def getDistanceMeters(lat1: Double, lon1: Double, lat2: Double, lon2: Double) = {
val theta = lon1 - lon2
Math.abs(
Math.round(
rad2deg(Math.acos(dist)) * 60 * 1.1515 * 1.609344 * 1000)
)
}
``````

Here's my implementation in Elixir

``````defmodule Geo do
@feet_per_sm 5280

@d2r :math.pi / 180

def deg_to_rad(deg), do: deg * @d2r

def great_circle_distance(p1, p2, :km), do: haversine(p1, p2) * @earth_radius_km
def great_circle_distance(p1, p2, :sm), do: haversine(p1, p2) * @earth_radius_sm
def great_circle_distance(p1, p2, :nm), do: haversine(p1, p2) * @earth_radius_nm
def great_circle_distance(p1, p2, :m), do: great_circle_distance(p1, p2, :km) * 1000
def great_circle_distance(p1, p2, :ft), do: great_circle_distance(p1, p2, :sm) * @feet_per_sm

@doc """
Calculate the [Haversine](https://en.wikipedia.org/wiki/Haversine_formula)
distance between two coordinates. Result is in radians. This result can be
multiplied by the sphere's radius in any unit to get the distance in that unit.
For example, multiple the result of this function by the Earth's radius in
kilometres and you get the distance between the two given points in kilometres.
"""
def haversine({lat1, lon1}, {lat2, lon2}) do
dlat = deg_to_rad(lat2 - lat1)
dlon = deg_to_rad(lon2 - lon1)

a = :math.pow(:math.sin(dlat / 2), 2) +
:math.pow(:math.sin(dlon / 2), 2) *

2 * :math.atan2(:math.sqrt(a), :math.sqrt(1 - a))
end
end
``````

In Python, you can use the geopy library to compute the geodesic distance using the WGS84 ellipsoid:

``````from geopy.distance import geodesic
newport_ri = (41.49008, -71.312796)
cleveland_oh = (41.499498, -81.695391)
print(geodesic(newport_ri, cleveland_oh).km)
``````

Dart Version

Haversine Algorithm.

``````import 'dart:math';

class GeoUtils {

static double _degreesToRadians(degrees) {
return degrees * pi / 180;
}

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

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

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

I think a version of the algorithm in R is still missing:

``````gpsdistance<-function(lat1,lon1,lat2,lon2){

# internal function to change deg to rad

degreesToRadians<- function (degrees) {
return (degrees * pi / 180)
}

R<-6371e3  #radius of Earth in meters

phi1<-degreesToRadians(lat1) # latitude 1
phi2<-degreesToRadians(lat2) # latitude 2
lambda1<-degreesToRadians(lon1) # longitude 1
lambda2<-degreesToRadians(lon2) # longitude 2

delta_phi<-phi1-phi2 # latitude-distance
delta_lambda<-lambda1-lambda2 # longitude-distance

a<-sin(delta_phi/2)*sin(delta_phi/2)+
cos(phi1)*cos(phi2)*sin(delta_lambda/2)*
sin(delta_lambda/2)

cc<-2*atan2(sqrt(a),sqrt(1-a))

distance<- R * cc

return(distance)  # in meters
}
``````

For java

``````public static double degreesToRadians(double degrees) {
return degrees * Math.PI / 180;
}

public static double distanceInKmBetweenEarthCoordinates(Location location1, Location location2) {
double earthRadiusKm = 6371;

double dLat = degreesToRadians(location2.getLatitude()-location1.getLatitude());
double dLon = degreesToRadians(location2.getLongitude()-location1.getLongitude());

double lat1 = degreesToRadians(location1.getLatitude());
double lat2 = degreesToRadians(location2.getLatitude());

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));
return earthRadiusKm * c;
}
``````

For anyone searching for a Delphi/Pascal version:

``````function GreatCircleDistance(const Lat1, Long1, Lat2, Long2: Double): Double;
var
{"ft_lang_label":"__label__en","ft_lang_prob":0.68392426,"math_prob":0.9861473,"size":27168,"snap":"2022-40-2023-06","text_gpt3_token_len":8537,"char_repetition_ratio":0.14419821,"word_repetition_ratio":0.10801647,"special_character_ratio":0.3322659,"punctuation_ratio":0.17895752,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9988459,"pos_list":[0],"im_url_duplicate_count":[null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2022-10-06T10:54:21Z\",\"WARC-Record-ID\":\"<urn:uuid:767dc1c3-dc27-43b4-8483-7797738a293c>\",\"Content-Length\":\"463452\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:4e3140de-fdac-4410-8919-a2c19b7bb289>\",\"WARC-Concurrent-To\":\"<urn:uuid:3b096b1a-4101-4e62-8f8a-d10280ff156d>\",\"WARC-IP-Address\":\"151.101.1.69\",\"WARC-Target-URI\":\"https://stackoverflow.com/questions/365826/calculate-distance-between-2-gps-coordinates\",\"WARC-Payload-Digest\":\"sha1:KIGBXOMCIBE6RMP6ZQOKNHDFUKOH6YP6\",\"WARC-Block-Digest\":\"sha1:53RSO7XGD5VN6N7MU2B3KYLHBCYHR7G7\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2022/CC-MAIN-2022-40/CC-MAIN-2022-40_segments_1664030337803.86_warc_CC-MAIN-20221006092601-20221006122601-00104.warc.gz\"}"}