# Calculate distance between two latitude-longitude points? (Haversine formula)

962

How do I calculate the distance between two points specified by latitude and longitude?

For clarification, I'd like the distance in kilometers; the points use the WGS84 system and I'd like to understand the relative accuracies of the approaches available.

This question is tagged with algorithm math maps latitude-longitude haversine

1206

This link might be helpful to you, as it details the use of the Haversine formula to calculate the distance.

Excerpt:

This script [in Javascript] calculates great-circle distances between the two points – that is, the shortest distance over the earth’s surface – using the ‘Haversine’ formula.

function getDistanceFromLatLonInKm(lat1,lon1,lat2,lon2) {
var R = 6371; // Radius of the earth in km
var a =
Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.sin(dLon/2) * Math.sin(dLon/2)
;
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c; // Distance in km
return d;
}

return deg * (Math.PI/180)
}

405

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, pi

def distance(lat1, lon1, lat2, lon2):
p = pi/180
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)) #2*R*asin...

And for the sake of completeness: Haversine on wiki.

71

Here is a C# Implementation:

static class DistanceAlgorithm
{
const double PIx = 3.141592653589793;

/// <summary>
/// </summary>
/// <param name="x">Degrees</param>
{
return x * PIx / 180;
}

/// <summary>
/// Calculate the distance between two places.
/// </summary>
/// <param name="lon1"></param>
/// <param name="lat1"></param>
/// <param name="lon2"></param>
/// <param name="lat2"></param>
/// <returns></returns>
public static double DistanceBetweenPlaces(
double lon1,
double lat1,
double lon2,
double lat2)
{
double dlon = Radians(lon2 - lon1);
double dlat = Radians(lat2 - lat1);

double a = (Math.Sin(dlat / 2) * Math.Sin(dlat / 2)) + Math.Cos(Radians(lat1)) * Math.Cos(Radians(lat2)) * (Math.Sin(dlon / 2) * Math.Sin(dlon / 2));
double angle = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
}

}

68

Here is a java implementation of the Haversine formula.

public final static double AVERAGE_RADIUS_OF_EARTH_KM = 6371;
public int calculateDistanceInKilometer(double userLat, double userLng,
double venueLat, double venueLng) {

double latDistance = Math.toRadians(userLat - venueLat);
double lngDistance = Math.toRadians(userLng - venueLng);

double a = Math.sin(latDistance / 2) * Math.sin(latDistance / 2)
* Math.sin(lngDistance / 2) * Math.sin(lngDistance / 2);

double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

}

Note that here we are rounding the answer to the nearest km.

42

Thanks very much for all this. I used the following code in my Objective-C iPhone app:

const double PIx = 3.141592653589793;
const double RADIO = 6371; // Mean radius of Earth in Km

return val * PIx / 180;
}

-(double)kilometresBetweenPlace1:(CLLocationCoordinate2D) place1 andPlace2:(CLLocationCoordinate2D) place2 {

double dlon = convertToRadians(place2.longitude - place1.longitude);
double dlat = convertToRadians(place2.latitude - place1.latitude);

double a = ( pow(sin(dlat / 2), 2) + cos(convertToRadians(place1.latitude))) * cos(convertToRadians(place2.latitude)) * pow(sin(dlon / 2), 2);
double angle = 2 * asin(sqrt(a));

}

Latitude and Longitude are in decimal. I didn't use min() for the asin() call as the distances that I'm using are so small that they don't require it.

It gave incorrect answers until I passed in the values in Radians - now it's pretty much the same as the values obtained from Apple's Map app :-)

Extra update:

If you are using iOS4 or later then Apple provide some methods to do this so the same functionality would be achieved with:

-(double)kilometresBetweenPlace1:(CLLocationCoordinate2D) place1 andPlace2:(CLLocationCoordinate2D) place2 {

MKMapPoint  start, finish;

start = MKMapPointForCoordinate(place1);
finish = MKMapPointForCoordinate(place2);

return MKMetersBetweenMapPoints(start, finish) / 1000;
}

40

This is a simple PHP function that will give a very reasonable approximation (under +/-1% error margin).

<?php
function distance(\$lat1, \$lon1, \$lat2, \$lon2) {

\$pi80 = M_PI / 180;
\$lat1 *= \$pi80;
\$lon1 *= \$pi80;
\$lat2 *= \$pi80;
\$lon2 *= \$pi80;

\$r = 6372.797; // mean radius of Earth in km
\$dlat = \$lat2 - \$lat1;
\$dlon = \$lon2 - \$lon1;
\$a = sin(\$dlat / 2) * sin(\$dlat / 2) + cos(\$lat1) * cos(\$lat2) * sin(\$dlon / 2) * sin(\$dlon / 2);
\$c = 2 * atan2(sqrt(\$a), sqrt(1 - \$a));
\$km = \$r * \$c;

//echo '<br/>'.\$km;
return \$km;
}
?>

As said before; the earth is NOT a sphere. It is like an old, old baseball that Mark McGwire decided to practice with - it is full of dents and bumps. The simpler calculations (like this) treat it like a sphere.

Different methods may be more or less precise according to where you are on this irregular ovoid AND how far apart your points are (the closer they are the smaller the absolute error margin). The more precise your expectation, the more complex the math.

31

I post here my working example.

List all points in table having distance between a designated point (we use a random point - lat:45.20327, long:23.7806) less than 50 KM, with latitude & longitude, in MySQL (the table fields are coord_lat and coord_long):

List all having DISTANCE<50, in Kilometres (considered Earth radius 6371 KM):

FROM obiective
WHERE coord_lat<>''
AND coord_long<>''
HAVING distanta<50
ORDER BY distanta desc

The above example was tested in MySQL 5.0.95 and 5.5.16 (Linux).

30

In the other answers an implementation in is missing.

Calculating the distance between two point is quite straightforward with the distm function from the geosphere package:

distm(p1, p2, fun = distHaversine)

where:

p1 = longitude/latitude for point(s)
p2 = longitude/latitude for point(s)
# type of distance calculation
fun = distCosine / distHaversine / distVincentySphere / distVincentyEllipsoid

As the earth is not perfectly spherical, the Vincenty formula for ellipsoids is probably the best way to calculate distances. Thus in the geosphere package you use then:

distm(p1, p2, fun = distVincentyEllipsoid)

Off course you don't necessarily have to use geosphere package, you can also calculate the distance in base R with a function:

hav.dist <- function(long1, lat1, long2, lat2) {
R <- 6371
diff.long <- (long2 - long1)
diff.lat <- (lat2 - lat1)
a <- sin(diff.lat/2)^2 + cos(lat1) * cos(lat2) * sin(diff.long/2)^2
b <- 2 * asin(pmin(1, sqrt(a)))
d = R * b
return(d)
}

11

The haversine is definitely a good formula for probably most cases, other answers already include it so I am not going to take the space. But it is important to note that no matter what formula is used (yes not just one). Because of the huge range of accuracy possible as well as the computation time required. The choice of formula requires a bit more thought than a simple no brainer answer.

This posting from a person at nasa, is the best one I found at discussing the options

http://www.cs.nyu.edu/visual/home/proj/tiger/gisfaq.html

For example, if you are just sorting rows by distance in a 100 miles radius. The flat earth formula will be much faster than the haversine.

HalfPi = 1.5707963;
R = 3956; /* the radius gives you the measurement unit*/

u = a * a + b * b;
v = - 2 * a * b * cos(longdestrad - longoriginrad);
c = sqrt(abs(u + v));
return R * c;

Notice there is just one cosine and one square root. Vs 9 of them on the Haversine formula.

8

All the above answers assumes the earth is a sphere. However, a more accurate approximation would be that of an oblate spheroid.

def Distance(lat1, lons1, lat2, lons2):
R1=(((((a**2)*math.cos(lat1))**2)+(((b**2)*math.sin(lat1))**2))/((a*math.cos(lat1))**2+(b*math.sin(lat1))**2))**0.5 #radius of earth at lat1
x1=R*math.cos(lat1)*math.cos(lons1)
y1=R*math.cos(lat1)*math.sin(lons1)
z1=R*math.sin(lat1)

R1=(((((a**2)*math.cos(lat2))**2)+(((b**2)*math.sin(lat2))**2))/((a*math.cos(lat2))**2+(b*math.sin(lat2))**2))**0.5 #radius of earth at lat2
x2=R*math.cos(lat2)*math.cos(lons2)
y2=R*math.cos(lat2)*math.sin(lons2)
z2=R*math.sin(lat2)

return ((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)**0.5

8

Python implimentation Origin is the center of the contiguous United States.

from haversine import haversine
origin = (39.50, 98.35)
paris = (48.8567, 2.3508)
haversine(origin, paris, miles=True)

To get the answer in kilometers simply set miles=false.

8

There could be a simpler solution, and more correct: The perimeter of earth is 40,000Km at the equator, about 37,000 on Greenwich (or any longitude) cycle. Thus:

pythagoras = function (lat1, lon1, lat2, lon2) {
function sqr(x) {return x * x;}
function cosDeg(x) {return Math.cos(x * Math.PI / 180.0);}

var earthCyclePerimeter = 40000000.0 * cosDeg((lat1 + lat2) / 2.0);
var dx = (lon1 - lon2) * earthCyclePerimeter / 360.0;
var dy = 37000000.0 * (lat1 - lat2) / 360.0;

return Math.sqrt(sqr(dx) + sqr(dy));
};

I agree that it should be fine-tuned as, I myself said that it's an ellipsoid, so the radius to be multiplied by the cosine varies. But it's a bit more accurate. Compared with Google Maps and it did reduce the error significantly.

8

I don't like adding yet another answer, but the Google maps API v.3 has spherical geometry (and more). After converting your WGS84 to decimal degrees you can do this:

No word about how accurate Google's calculations are or even what model is used (though it does say "spherical" rather than "geoid". By the way, the "straight line" distance will obviously be different from the distance if one travels on the surface of the earth which is what everyone seems to be presuming.

7

You can use the build in CLLocationDistance to calculate this:

CLLocation *location1 = [[CLLocation alloc] initWithLatitude:latitude1 longitude:longitude1];
CLLocation *location2 = [[CLLocation alloc] initWithLatitude:latitude2 longitude:longitude2];
[self distanceInMetersFromLocation:location1 toLocation:location2]

- (int)distanceInMetersFromLocation:(CLLocation*)location1 toLocation:(CLLocation*)location2 {
CLLocationDistance distanceInMeters = [location1 distanceFromLocation:location2];
return distanceInMeters;
}

In your case if you want kilometers just divide by 1000.

6

As pointed out, an accurate calculation should take into account that the earth is not a perfect sphere. Here are some comparisons of the various algorithms offered here:

geoDistance(50,5,58,3)
Haversine: 899 km
Maymenn: 833 km
Keerthana: 897 km

geoDistance(50,5,-58,-3)
Haversine: 12030 km
Maymenn: 11135 km
Keerthana: 10310 km

geoDistance(.05,.005,.058,.003)
Haversine: 0.9169 km
Maymenn: 0.851723 km
Keerthana: 0.917964 km

geoDistance(.05,80,.058,80.3)
Haversine: 33.37 km
Maymenn: 33.34 km
Keerthana: 33.40767 km

Over small distances, Keerthana's algorithm does seem to coincide with that of Google Maps. Google Maps does not seem to follow any simple algorithm, suggesting that it may be the most accurate method here.

Anyway, here is a Javascript implementation of Keerthana's algorithm:

function geoDistance(lat1, lng1, lat2, lng2){
const a = 6378.137; // equitorial radius in km
const b = 6356.752; // polar radius in km

var sq = x => (x*x);
var sqr = x => Math.sqrt(x);
var cos = x => Math.cos(x);
var sin = x => Math.sin(x);
var radius = lat => sqr((sq(a*a*cos(lat))+sq(b*b*sin(lat)))/(sq(a*cos(lat))+sq(b*sin(lat))));

lat1 = lat1 * Math.PI / 180;
lng1 = lng1 * Math.PI / 180;
lat2 = lat2 * Math.PI / 180;
lng2 = lng2 * Math.PI / 180;

var x1 = R1*cos(lat1)*cos(lng1);
var y1 = R1*cos(lat1)*sin(lng1);
var z1 = R1*sin(lat1);

var x2 = R2*cos(lat2)*cos(lng2);
var y2 = R2*cos(lat2)*sin(lng2);
var z2 = R2*sin(lat2);

return sqr(sq(x1-x2)+sq(y1-y2)+sq(z1-z2));
}

6

Here is the SQL Implementation to calculate the distance in km,

SELECT UserId, ( 3959 * acos( cos( radians( your latitude here ) ) * cos( radians(latitude) ) *
sin( radians(latitude) ) ) ) AS distance FROM user HAVING
distance < 5  ORDER BY distance LIMIT 0 , 5;

For further details in the implementation by programming langugage, you can just go through the php script given here

5

Here is a typescript implementation of the Haversine formula

static getDistanceFromLatLonInKm(lat1: number, lon1: number, lat2: number, lon2: number): number {
var deg2Rad = deg => {
return deg * Math.PI / 180;
}

var r = 6371; // Radius of the earth in km
var dLat = deg2Rad(lat2 - lat1);
var dLon = deg2Rad(lon2 - lon1);
var a =
Math.sin(dLat / 2) * Math.sin(dLat / 2) +
Math.sin(dLon / 2) * Math.sin(dLon / 2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
var d = r * c; // Distance in km
return d;
}

4

This script [in PHP] calculates distances between the two points.

public static function getDistanceOfTwoPoints(\$source, \$dest, \$unit='K') {
\$lat1 = \$source[0];
\$lon1 = \$source[1];
\$lat2 = \$dest[0];
\$lon2 = \$dest[1];

\$theta = \$lon1 - \$lon2;
\$dist = acos(\$dist);
\$miles = \$dist * 60 * 1.1515;
\$unit = strtoupper(\$unit);

if (\$unit == "K") {
return (\$miles * 1.609344);
}
else if (\$unit == "M")
{
return (\$miles * 1.609344 * 1000);
}
else if (\$unit == "N") {
return (\$miles * 0.8684);
}
else {
return \$miles;
}
}

4

Java implementation in according Haversine formula

double calculateDistance(double latPoint1, double lngPoint1,
double latPoint2, double lngPoint2) {
if(latPoint1 == latPoint2 && lngPoint1 == lngPoint2) {
return 0d;
}

final double EARTH_RADIUS = 6371.0; //km value;

double distance = Math.pow(Math.sin((latPoint2 - latPoint1) / 2.0), 2)
+ Math.cos(latPoint1) * Math.cos(latPoint2)
* Math.pow(Math.sin((lngPoint2 - lngPoint1) / 2.0), 2);
distance = 2.0 * EARTH_RADIUS * Math.asin(Math.sqrt(distance));

return distance; //km value
}

3

Here's the accepted answer implementation ported to Java in case anyone needs it.

package com.project529.garage.util;

/**
*/
private static double EARTH_RADIUS = 6371;

/**
* Returns the distance between two sets of latitudes and longitudes in meters.
* <p/>
* Based from the following JavaScript SO answer:
* http://stackoverflow.com/questions/27928/calculate-distance-between-two-latitude-longitude-points-haversine-formula,
* which is based on https://en.wikipedia.org/wiki/Haversine_formula (error rate: ~0.55%).
*/
public double getDistanceBetween(double lat1, double lon1, double lat2, double lon2) {
double dLat = toRadians(lat2 - lat1);
double dLon = toRadians(lon2 - lon1);

double a = Math.sin(dLat / 2) * Math.sin(dLat / 2) +
Math.sin(dLon / 2) * Math.sin(dLon / 2);
double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
double d = EARTH_RADIUS * c;

return d;
}

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

3

To calculate the distance between two points on a sphere you need to do the Great Circle calculation.

There are a number of C/C++ libraries to help with map projection at MapTools if you need to reproject your distances to a flat surface. To do this you will need the projection string of the various coordinate systems.

You may also find MapWindow a useful tool to visualise the points. Also as its open source its a useful guide to how to use the proj.dll library, which appears to be the core open source projection library.

3

Here is my java implementation for calculation distance via decimal degrees after some search. I used mean radius of world (from wikipedia) in km. If you want result miles then use world radius in miles.

public static double distanceLatLong2(double lat1, double lng1, double lat2, double lng2)
{
double earthRadius = 6371.0d; // KM: use mile here if you want mile result

double dLat = toRadian(lat2 - lat1);
double dLng = toRadian(lng2 - lng1);

double a = Math.pow(Math.sin(dLat/2), 2)  +
Math.pow(Math.sin(dLng/2), 2);

double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));

return earthRadius * c; // returns result kilometers
}

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

3

For those looking for an Excel formula based on WGS-84 & GRS-80 standards:

Source

2

Here's another converted to Ruby code:

include Math
#Note: from/to = [lat, long]

def get_distance_in_km(from, to)
radians = lambda { |deg| deg * Math.PI / 180 }

cosines_product = Math.sin(dLat/2) * Math.sin(dLat/2) + Math.cos(radians[from[0]]) * Math.cos(radians[to[1]]) * Math.sin(dLon/2) * Math.sin(dLon/2)

c = 2 * Math.atan2(Math.sqrt(cosines_product), Math.sqrt(1-cosines_product))
return radius * c # Distance in kilometer
end

2

function getDistanceFromLatLonInKm(position1, position2) {
"use strict";
var deg2rad = function (deg) { return deg * (Math.PI / 180); },
R = 6371,
a = Math.sin(dLat / 2) * Math.sin(dLat / 2)
* Math.sin(dLng / 2) * Math.sin(dLng / 2),
c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
return R * c;
}

console.log(getDistanceFromLatLonInKm(
{lat: 48.7931459, lng: 1.9483572},
{lat: 48.827167, lng: 2.2459745}
));

2

Here is the implementation VB.NET, this implementation will give you the result in KM or Miles based on an Enum value you pass.

Public Enum DistanceType
Miles
KiloMeters
End Enum

Public Structure Position
Public Latitude As Double
Public Longitude As Double
End Structure

Public Class Haversine

Public Function Distance(Pos1 As Position,
Pos2 As Position,
DistType As DistanceType) As Double

Dim R As Double = If((DistType = DistanceType.Miles), 3960, 6371)

Dim dLat As Double = Me.toRadian(Pos2.Latitude - Pos1.Latitude)

Dim dLon As Double = Me.toRadian(Pos2.Longitude - Pos1.Longitude)

Dim a As Double = Math.Sin(dLat / 2) * Math.Sin(dLat / 2) + Math.Cos(Me.toRadian(Pos1.Latitude)) * Math.Cos(Me.toRadian(Pos2.Latitude)) * Math.Sin(dLon / 2) * Math.Sin(dLon / 2)

Dim c As Double = 2 * Math.Asin(Math.Min(1, Math.Sqrt(a)))

Dim result As Double = R * c

Return result

End Function

Private Function toRadian(val As Double) As Double

Return (Math.PI / 180) * val

End Function

End Class

2

here is an example in postgres sql (in km, for miles version, replace 1.609344 by 0.8684 version)

CREATE OR REPLACE FUNCTION public.geodistance(alat float, alng float, blat

float, blng  float)
RETURNS float AS
\$BODY\$
DECLARE
v_distance float;
BEGIN

v_distance = asin( sqrt(
+ (
)
)
) * cast('7926.3352' as float) * cast('1.609344' as float) ;

RETURN v_distance;
END
\$BODY\$
language plpgsql VOLATILE SECURITY DEFINER;
alter function geodistance(alat float, alng float, blat float, blng float)
owner to postgres;

2

function getDistanceFromLatLonInKm(lat1,lon1,lat2,lon2,units) {
var R = 6371; // Radius of the earth in km
var a =
Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.sin(dLon/2) * Math.sin(dLon/2)
;
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c;
var miles = d / 1.609344;

if ( units == 'km' ) {
return d;
} else {
return miles;
}}

Chuck's solution, valid for miles also.

2

In Mysql use the following function pass the parameters as using POINT(LONG,LAT)

CREATE FUNCTION `distance`(a POINT, b POINT)
RETURNS double
DETERMINISTIC
BEGIN

RETURN

GLength( LineString(( PointFromWKB(a)), (PointFromWKB(b)))) * 100000; -- To Make the distance in meters

END;

2

I condensed the computation down by simplifying the formula.

Here it is in Ruby:

include Math
radians = lambda { |deg| deg * PI / 180 }

# from/to = { :lat => (latitude_in_degrees), :lng => (longitude_in_degrees) }
def haversine_distance(from, to)
cosines_product = cos(to[:lat]) * cos(from[:lat]) * cos(from[:lng] - to[:lng])
sines_product = sin(to[:lat]) * sin(from[:lat])
return earth_radius_mi * acos(cosines_product + sines_product)
end

2

As this is the most popular discussion of the topic I'll add my experience from late 2019-early 2020 here. To add to the existing answers - my focus was to find an accurate AND fast (i.e. vectorized) solution.

Let's start with what is mostly used by answers here - the Haversine approach. It is trivial to vectorize, see example in python below:

def haversine(lat1, lon1, lat2, lon2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)

All args must be of equal length.
Distances are in meters.

Ref:
https://stackoverflow.com/questions/29545704/fast-haversine-approximation-python-pandas
"""
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])

dlon = lon2 - lon1
dlat = lat2 - lat1

a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2

c = 2 * np.arcsin(np.sqrt(a))

# initial azimuth in degrees
y = np.sin(lon2-lon1) * np.cos(lat2)
x = np.cos(lat1)*np.sin(lat2) - np.sin(lat1)*np.cos(lat2)*np.cos(dlon)
azi1 = np.arctan2(y, x)*180./math.pi

return {'s12':s12, 'azi1': azi1}

Accuracy-wise, it is least accurate. Wikipedia states 0.5% of relative deviation on average without any sources. My experiments show less of a deviation. Below is the comparison ran on 100,000 random points vs my library, which should be accurate to millimeter levels:

np.random.seed(42)
lats1 = np.random.uniform(-90,90,100000)
lons1 = np.random.uniform(-180,180,100000)
lats2 = np.random.uniform(-90,90,100000)
lons2 = np.random.uniform(-180,180,100000)
r1 = inverse(lats1, lons1, lats2, lons2)
r2 = haversine(lats1, lons1, lats2, lons2)
print("Max absolute error: {:4.2f}m".format(np.max(r1['s12']-r2['s12'])))
print("Mean absolute error: {:4.2f}m".format(np.mean(r1['s12']-r2['s12'])))
print("Max relative error: {:4.2f}%".format(np.max((r2['s12']/r1['s12']-1)*100)))
print("Mean relative error: {:4.2f}%".format(np.mean((r2['s12']/r1['s12']-1)*100)))

Output:

Max absolute error: 26671.47m
Mean absolute error: -2499.84m
Max relative error: 0.55%
Mean relative error: -0.02%

So on average 2.5km deviation on 100,000 random pairs of coordinates, which may be good for majority of cases.

Next option is Vincenty's formulae which is accurate up to millimeters, depending on convergence criteria and can be vectorized as well. It does have the issue with convergence near antipodal points. You can make it converge at those points by relaxing convergence criteria, but accuracy drops to 0.25% and more. Outside of antipodal points Vincenty will provide results close to Geographiclib within relative error of less than 1.e-6 on average.

Geographiclib, mentioned here, is really the current golden standard. It has several implementations and fairly fast, especially if you are using C++ version.

Now, if you are planning to use Python for anything above 10k points I'd suggest to consider my vectorized implementation. I created a geovectorslib library with vectorized Vincenty routine for my own needs, which uses Geographiclib as fallback for near antipodal points. Below is the comparison vs Geographiclib for 100k points. As you can see it provides up to 20x improvement for inverse and 100x for direct methods for 100k points and the gap will grow with number of points. Accuracy-wise it will be within 1.e-5 rtol of Georgraphiclib.

Direct method for 100,000 points
94.9 ms ± 25 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
9.79 s ± 1.4 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

Inverse method for 100,000 points
1.5 s ± 504 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
24.2 s ± 3.91 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

2

I made a custom function in R to calculate haversine distance(km) between two spatial points using functions available in R base package.

custom_hav_dist <- function(lat1, lon1, lat2, lon2) {
R <- 6371

distance_in_km <- 6371*acos((cos(lat_1)*cos(lat_2))+
(sin(lat_1)*sin(lat_2)*cos(diff_long)))
rm(lat1, lon1, lat2, lon2)
return(distance_in_km)
}

Sample output

custom_hav_dist(50.31,19.08,54.14,19.39)
[1] 426.3987

PS: To calculate distances in miles, substitute R in function (6371) with 3958.756 (and for nautical miles, use 3440.065).

1

FSharp version, using miles:

let radialDistanceHaversine location1 location2 : float =
let degreeToRadian degrees = degrees * System.Math.PI / 180.0
let deltaLat = location2.Latitude - location1.Latitude |> degreeToRadian
let deltaLong = location2.Longitude - location1.Longitude |> degreeToRadian
let a =
(deltaLat / 2.0 |> sin) ** 2.0
+ (location1.Latitude |> degreeToRadian |> cos)
* (location2.Latitude |> degreeToRadian |> cos)
* (deltaLong / 2.0 |> sin) ** 2.0
atan2 (a |> sqrt) (1.0 - a |> sqrt)
* 2.0

1

One of the main challenges to calculating distances - especially large ones - is accounting for the curvature of the Earth. If only the Earth were flat, calculating the distance between two points would be as simple as for that of a straight line! The Haversine formula includes a constant (it's the R variable below) that represents the radius of the Earth. Depending on whether you are measuring in miles or kilometers, it would equal 3956 mi or 6367 km respectively.

The basic formula is:

dlon = lon2 - lon1
dlat = lat2 - lat1
a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
c = 2 * atan2( sqrt(a), sqrt(1-a) )
distance = R * c (where R is the radius of the Earth)

R = 6367 km OR 3956 mi
lat1, lon1: The Latitude and Longitude of point 1 (in decimal degrees)
lat2, lon2: The Latitude and Longitude of point 2 (in decimal degrees)
unit: The unit of measurement in which to calculate the results where:
'M' is statute miles (default)
'K' is kilometers
'N' is nautical miles

Sample

function distance(lat1, lon1, lat2, lon2, unit) {
try {
var radlat1 = Math.PI * lat1 / 180
var radlat2 = Math.PI * lat2 / 180
var theta = lon1 - lon2
var radtheta = Math.PI * theta / 180
dist = Math.acos(dist)
dist = dist * 180 / Math.PI
dist = dist * 60 * 1.1515
if (unit == "K") {
dist = dist * 1.609344
}
if (unit == "N") {
dist = dist * 0.8684
}
return dist
} catch (err) {
console.log(err);
}
}

1

there is a good example in here to calculate distance with PHP http://www.geodatasource.com/developers/php :

function distance(\$lat1, \$lon1, \$lat2, \$lon2, \$unit) {

\$theta = \$lon1 - \$lon2;
\$dist = acos(\$dist);
\$miles = \$dist * 60 * 1.1515;
\$unit = strtoupper(\$unit);

if (\$unit == "K") {
return (\$miles * 1.609344);
} else if (\$unit == "N") {
return (\$miles * 0.8684);
} else {
return \$miles;
}
}

1

Had an issue with math.deg in LUA... if anyone knows a fix please clean up this code!

In the meantime here's an implementation of the Haversine in LUA (use this with Redis!)

function calcDist(lat1, lon1, lat2, lon2)
lat1= lat1*0.0174532925
lat2= lat2*0.0174532925
lon1= lon1*0.0174532925
lon2= lon2*0.0174532925

dlon = lon2-lon1
dlat = lat2-lat1

a = math.pow(math.sin(dlat/2),2) + math.cos(lat1) * math.cos(lat2) * math.pow(math.sin(dlon/2),2)
c = 2 * math.asin(math.sqrt(a))
dist = 6371 * c      -- multiply by 0.621371 to convert to miles
return dist
end

cheers!

1

Dart lang:

import 'dart:math' show cos, sqrt, asin;

double calculateDistance(LatLng l1, LatLng l2) {
const p = 0.017453292519943295;
final a = 0.5 -
cos((l2.latitude - l1.latitude) * p) / 2 +
cos(l1.latitude * p) *
cos(l2.latitude * p) *
(1 - cos((l2.longitude - l1.longitude) * p)) /
2;
return 12742 * asin(sqrt(a));
}

0

The functions needed for an accurate calculation of distance between lat-long points are complex, and the pitfalls are many. I would not recomend haversine or other spherical solutions due to the big inaccuracies (the earth is not a perfect sphere). The vincenty formula is better, but will in some cases throw errors, even when coded correctly.

Instead of coding the functions yourself I suggest using geopy which have implemented the very accurate geographiclib for distance calculations (paper from author).

#pip install geopy
from geopy.distance import geodesic
NY = [40.71278,-74.00594]
Beijing = [39.90421,116.40739]
print("WGS84: ",geodesic(NY, Beijing).km) #WGS84 is Standard
print("Intl24: ",geodesic(NY, Beijing, ellipsoid='Intl 1924').km) #geopy includes different ellipsoids
print("Custom ellipsoid: ",geodesic(NY, Beijing, ellipsoid=(6377., 6356., 1 / 297.)).km) #custom ellipsoid

#supported ellipsoids:
#model             major (km)   minor (km)     flattening
#'WGS-84':        (6378.137,    6356.7523142,  1 / 298.257223563)
#'GRS-80':        (6378.137,    6356.7523141,  1 / 298.257222101)
#'Airy (1830)':   (6377.563396, 6356.256909,   1 / 299.3249646)
#'Intl 1924':     (6378.388,    6356.911946,   1 / 297.0)
#'Clarke (1880)': (6378.249145, 6356.51486955, 1 / 293.465)
#'GRS-67':        (6378.1600,   6356.774719,   1 / 298.25)

The only drawback with this library is that it doesn't support vectorized calculations. For vectorized calculations you can use the new geovectorslib.

#pip install geovectorslib
from geovectorslib import inverse
print(inverse(lats1,lons1,lats2,lons2)['s12'])

lats and lons are numpy arrays. Geovectorslib is very accurate and extremly fast! I haven't found a solution for changing ellipsoids though. The WGS84 ellipsoid is used as standard, which is the best choice for most uses.

0

Here's a simple javascript function that may be useful from this link.. somehow related but we're using google earth javascript plugin instead of maps

function getApproximateDistanceUnits(point1, point2) {

var xs = 0;
var ys = 0;

xs = point2.getX() - point1.getX();
xs = xs * xs;

ys = point2.getY() - point1.getY();
ys = ys * ys;

return Math.sqrt(xs + ys);
}

The units tho are not in distance but in terms of a ratio relative to your coordinates. There are other computations related you can substitute for the getApproximateDistanceUnits function link here

Then I use this function to see if a latitude longitude is within the radius

if (point1 && point2) {
} else {
return 0;
}
}

point may be defined as

\$\$.getPoint = function(lati, longi) {
var location = {
x: 0,
y: 0,
getX: function() { return location.x; },
getY: function() { return location.y; }
};
location.x = lati;
location.y = longi;

return location;
};

then you can do your thing to see if a point is within a region with a radius say:

//put it on the map if within the range of a specified radi assuming 100,000,000 units
var iconpoint = Map.getPoint(pp.latitude, pp.longitude);
var centerpoint = Map.getPoint(Settings.CenterLatitude, Settings.CenterLongitude);

//approx ~200 units to show only half of the globe from the default center radius
}
else {
otherSidePlacemarks.push({
latitude: pp.latitude,
longitude: pp.longitude,
name: pp.name
});

}

0

If you want the driving distance/route (posting it here because this is the first result for the distance between two points on google but for most people the driving distance is more useful), you can use Google Maps Distance Matrix Service:

getDrivingDistanceBetweenTwoLatLong(origin, destination) {

return new Observable(subscriber => {
service.getDistanceMatrix(
{
travelMode: 'DRIVING'
}, (response, status) => {
console.log('Error:', status);
subscriber.error({error: status, status: status});
} else {
console.log(response);
try {
let valueInMeters = response.rows[0].elements[0].distance.value;
let valueInKms = valueInMeters / 1000;
subscriber.next(valueInKms);
subscriber.complete();
}
catch(error) {
subscriber.error({error: error, status: status});
}
}
});
});
}

0

If you are using python; pip install geopy

from geopy.distance import geodesic

origin = (30.172705, 31.526725)  # (latitude, longitude) don't confuse
destination = (30.288281, 31.732326)

print(geodesic(origin, destination).meters)  # 23576.805481751613
print(geodesic(origin, destination).kilometers)  # 23.576805481751613
print(geodesic(origin, destination).miles)  # 14.64994773134371

0

//JAVA
public Double getDistanceBetweenTwoPoints(Double latitude1, Double longitude1, Double latitude2, Double longitude2) {

double dLat = getRad(latitude2 - latitude1);
double dLong = getRad(longitude2 - longitude1);

double a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.cos(getRad(latitude1)) * Math.cos(getRad(latitude2)) * Math.sin(dLong / 2) * Math.sin(dLong / 2);
double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
return (RADIUS_EARTH * c) * 1000;
}

return x * Math.PI / 180;
}

0

I've created this small Javascript LatLng object, might be useful for somebody.

var latLng1 = new LatLng(5, 3);
var latLng2 = new LatLng(6, 7);
var distance = latLng1.distanceTo(latLng2);

Code:

/**
* latLng point
* @param {Number} lat
* @param {Number} lng
* @returns {LatLng}
* @constructor
*/
function LatLng(lat,lng) {
this.lat = parseFloat(lat);
this.lng = parseFloat(lng);

this.__cache = {};
}

LatLng.prototype = {
toString: function() {
return [this.lat, this.lng].join(",");
},

/**
* calculate distance in km to another latLng, with caching
* @param {LatLng} latLng
* @returns {Number} distance in km
*/
distanceTo: function(latLng) {
var cacheKey = latLng.toString();
if(cacheKey in this.__cache) {
return this.__cache[cacheKey];
}

// the fastest way to calculate the distance, according to this jsperf test;
// http://stackoverflow.com/questions/27928
var deg2rad = 0.017453292519943295; // === Math.PI / 180
var lat1 = this.lat * deg2rad;
var lng1 = this.lng * deg2rad;
var lat2 = latLng.lat * deg2rad;
var lng2 = latLng.lng * deg2rad;
var a = (
(1 - Math.cos(lat2 - lat1)) +
(1 - Math.cos(lng2 - lng1)) * Math.cos(lat1) * Math.cos(lat2)
) / 2;
var distance = 12742 * Math.asin(Math.sqrt(a)); // Diameter of the earth in km (2 * 6371)

// cache the distance
this.__cache[cacheKey] = distance;

return distance;
}
};

0

You can calculate it by using Haversine formula which is:

a = sin²(?f/2) + cos f1 · cos f2 · sin²(??/2)
c = 2 · atan2( va, v(1-a) )
d = R · c

An example to calculate distance between two points is given below

Suppose i have to calculate distance between New Delhi to London, so how can i use this formula :

New delhi co-ordinates= 28.7041° N, 77.1025° E
London co-ordinates= 51.5074° N, 0.1278° W

var R = 6371e3; // metres