Reputation: 15357
Calculate distance between two latitude-longitude points? (Haversine formula)
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.
Upvotes: 1172
Views: 1232254
Answers (30)
Reputation: 36
This is a direct python implementation of the Haversine formula. This implementation returns the distance between the two points in meters.
import numpy as np
def distance_from_lat_long(lat1, lon1, lat2, lon2):
R = 6371000 # earth radius in meters
# convert lat and lon to radians
lat1 = np.deg2rad(lat1)
lat2 = np.deg2rad(lat2)
lon1 = np.deg2rad(lon1)
lon2 = np.deg2rad(lon2)
# haversine formula
return (2 * R * np.arcsin(
np.sqrt((1 - np.cos(lat2 - lat1) + np.cos(lat1) * np.cos(lat2) *
(1 - np.cos(lon2 - lon1))) / 2)))
Upvotes: 1
Reputation: 83215
In the other answers an implementation in r 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
p <- pi/180
diff.long <- (long2 - long1) * p
diff.lat <- (lat2 - lat1) * p
a <- sin(diff.lat/2)^2 + cos(lat1 * p) * cos(lat2 * p) * sin(diff.long/2)^2
b <- 2 * asin(pmin(1, sqrt(a)))
d = R * b
return(d)
}
Upvotes: 32
Reputation: 2072
To calculate the distance between two points specified by latitude and longitude, you can use the Haversine formula or Vincenty's formulae. Both formulas take into account the curvature of the Earth and provide reasonably accurate results. Here's a brief explanation of each method:
Haversine Formula: The Haversine formula calculates the great-circle distance between two points on a sphere (in this case, the Earth) assuming a perfect sphere. It is simpler and faster than Vincenty's formulae but may have slightly lower accuracy, especially for longer distances.
The Haversine formula is as follows:
a = sin²(Δlat/2) + cos(lat1) * cos(lat2) * sin²(Δlon/2)
c = 2 * atan2(√a, √(1-a))
d = R * c
Where:
Δlat is the difference in latitude between the two points.
Δlon is the difference in longitude between the two points.
lat1 and lat2 are the latitudes of the two points.
R is the radius of the Earth (mean radius: 6,371 km).
You can find direct implementation at Haversine Implementations
Vincenty's Formulae: Vincenty's formulae are more complex but provide higher accuracy for calculating the distance between two points on an ellipsoidal model of the Earth. They take into account the Earth's shape and provide accurate results for any pair of points on the globe.
There are two variations of Vincenty's formulae: the direct formula, which calculates the distance between points, and the inverse formula, which calculates the initial bearing, final bearing, and distance between points.
To use Vincenty's formulae, you'll need to implement the specific equations in your chosen programming language.
Direct Formula: The direct formula calculates the destination point (latitude and longitude) given a starting point, initial bearing, and distance. It is useful when you know the starting point, the desired distance, and the initial bearing towards the destination.
The formula is as follows:
α1 = atan2((1 - f) * tan(lat1), cos(α0))
sinσ1 = (cos(U2) * sin(α1))^2 + (cos(U1) * sin(U2) - sin(U1) * cos(U2) * cos(α1))^2
sinσ = sqrt(sinσ1)
cosσ = sin(U1) * sin(U2) + cos(U1) * cos(U2) * cos(α1)
σ = atan2(sinσ, cosσ)
sinα = cos(U1) * cos(U2) * sin(α1) / sinσ
cos2αm = 1 - sinα^2
C = (f / 16) * cos2αm * (4 + f * (4 - 3 * cos2αm))
λ = L + (1 - C) * f * sinα * (σ + C * sinσ * (cos(2 * σ) + C * cosσ * (-1 + 2 * cos(2 * σ)^2)))
Δ = L2 - λ
L = λ
Where:
lat1, lon1: Starting point latitude and longitude in radians.
α0: Initial bearing in radians.
U1 = atan((1 - f) * tan(lat1))
U2 = atan((1 - f) * tan(lat2)), where lat2 is calculated iteratively using the formula above.
f = (a - b) / a, where a and b are the equatorial and polar radii of the Earth, respectively.
Inverse Formula: The inverse formula calculates the distance, initial bearing, and final bearing between two points on the Earth's surface. It is useful when you know the latitude and longitude of both points.
The formula is as follows:
L = lon2 - lon1
λ = L
λʹ = 2π + atan2(y, x)
σ = atan2(yʹ, xʹ)
C = (f / 16) * cos^2αm * (4 + f * (4 - 3 * cos^2αm))
λ = L + (1 - C) * f * sinα * (σ + C * sinσ * (cos(2 * σ) + C * cosσ * (-1 + 2 * cos(2 * σ)^2)))
Where:
lat1, lon1: Latitude and longitude of the first point in radians.
lat2, lon2: Latitude and longitude of the second point in radians.
L = lon2 - lon1
U1 = atan((1 - f) * tan(lat1))
U2 = atan((1 - f) * tan(lat2))
sinα = cos(U2) * sin(L)
cosα = sqrt(1 - sinα^2)
cos^2αm = cosα^2 * cosσ^2
x = σ - sinα * sinσ
y = λʹ - sinα * sinσ
xʹ = cos(U1) * sin(U2) - sin(U1) * cos(U2) * cos(L)
yʹ = cos(U2) * sin(L)
Upvotes: 1
Reputation: 222441
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) {
const r = 6371; // km
const p = Math.PI / 180;
const a = 0.5 - Math.cos((lat2 - lat1) * p) / 2
+ Math.cos(lat1 * p) * Math.cos(lat2 * p) *
(1 - Math.cos((lon2 - lon1) * p)) / 2;
return 2 * r * Math.asin(Math.sqrt(a));
}
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):
r = 6371 # km
p = pi / 180
a = 0.5 - cos((lat2-lat1)*p)/2 + cos(lat1*p) * cos(lat2*p) * (1-cos((lon2-lon1)*p))/2
return 2 * r * asin(sqrt(a))
And for the sake of completeness: Haversine on Wikipedia.
Upvotes: 517
Reputation: 700
Java implementation in according to 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;
//converting to radians
latPoint1 = Math.toRadians(latPoint1);
lngPoint1 = Math.toRadians(lngPoint1);
latPoint2 = Math.toRadians(latPoint2);
lngPoint2 = Math.toRadians(lngPoint2);
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
}
Upvotes: 7
Reputation: 898
There is some errors in the code provided, I've fixed it below.
All the above answers assumes the earth is a sphere. However, a more accurate approximation would be that of an oblate spheroid.
a= 6378.137#equitorial radius in km
b= 6356.752#polar radius in km
def Distance(lat1, lons1, lat2, lons2):
lat1=math.radians(lat1)
lons1=math.radians(lons1)
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=R1*math.cos(lat1)*math.cos(lons1)
y1=R1*math.cos(lat1)*math.sin(lons1)
z1=R1*math.sin(lat1)
lat2=math.radians(lat2)
lons2=math.radians(lons2)
R2=(((((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=R2*math.cos(lat2)*math.cos(lons2)
y2=R2*math.cos(lat2)*math.sin(lons2)
z2=R2*math.sin(lat2)
return ((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)**0.5
Upvotes: 10
Reputation: 3572
Here's a Scala implementation:
def calculateHaversineDistance(lat1: Double, lon1: Double, lat2: Double, lon2: Double): Double = {
val long2 = lon2 * math.Pi / 180
val lat2 = lat2 * math.Pi / 180
val long1 = lon1 * math.Pi / 180
val lat1 = lat1 * math.Pi / 180
val dlon = long2 - long1
val dlat = lat2 - lat1
val a = math.pow(math.sin(dlat / 2), 2) + math.cos(lat1) * math.cos(lat2) * math.pow(math.sin(dlon / 2), 2)
val c = 2 * math.atan2(Math.sqrt(a), math.sqrt(1 - a))
val haversineDistance = 3961 * c // 3961 = radius of earth in miles
haversineDistance
}
Upvotes: 0
Reputation: 2428
You could use a module like geolib too:
How to install:
$ npm install geolib
How to use:
import { getDistance } from 'geolib'
const distance = getDistance(
{ latitude: 51.5103, longitude: 7.49347 },
{ latitude: "51° 31' N", longitude: "7° 28' E" }
)
console.log(distance)
Documentation: https://www.npmjs.com/package/geolib
Upvotes: 2
Reputation: 3897
pip install haversine
Python implementation
Origin is the center of the contiguous United States.
from haversine import haversine, Unit
origin = (39.50, 98.35)
paris = (48.8567, 2.3508)
haversine(origin, paris, unit=Unit.MILES)
To get the answer in kilometers simply set unit=Unit.KILOMETERS (that's the default).
Upvotes: 9
Reputation: 2402
function getDistanceFromLatLonInKm(position1, position2) {
"use strict";
var deg2rad = function (deg) { return deg * (Math.PI / 180); },
R = 6371,
dLat = deg2rad(position2.lat - position1.lat),
dLng = deg2rad(position2.lng - position1.lng),
a = Math.sin(dLat / 2) * Math.sin(dLat / 2)
+ Math.cos(deg2rad(position1.lat))
* Math.cos(deg2rad(position2.lat))
* 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}
));
Upvotes: 3
Reputation: 20760
Here is the Erlang implementation
lat_lng({Lat1, Lon1}=_Point1, {Lat2, Lon2}=_Point2) ->
P = math:pi() / 180,
R = 6371, % Radius of Earth in KM
A = 0.5 - math:cos((Lat2 - Lat1) * P) / 2 +
math:cos(Lat1 * P) * math:cos(Lat2 * P) * (1 - math:cos((Lon2 - Lon1) * P))/2,
R * 2 * math:asin(math:sqrt(A)).
Upvotes: 1
Reputation: 83
function distance($lat1, $lon1, $lat2, $lon2) {
$pi80 = M_PI / 180;
$lat1 *= $pi80; $lon1 *= $pi80; $lat2 *= $pi80; $lon2 *= $pi80;
$dlat = $lat2 - $lat1;
$dlon = $lon2 - $lon1;
$a = sin($dlat / 2) * sin($dlat / 2) + cos($lat1) * cos($lat2) * sin($dlon / 2) * sin($dlon / 2);
$km = 6372.797 * 2 * atan2(sqrt($a), sqrt(1 - $a));
return $km;
}
Upvotes: 0
Reputation: 6054
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 => {
let service = new google.maps.DistanceMatrixService();
service.getDistanceMatrix(
{
origins: [new google.maps.LatLng(origin.lat, origin.long)],
destinations: [new google.maps.LatLng(destination.lat, destination.long)],
travelMode: 'DRIVING'
}, (response, status) => {
if (status !== google.maps.DistanceMatrixStatus.OK) {
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});
}
}
});
});
}
Upvotes: 2
Reputation: 4428
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
var dist = Math.sin(radlat1) * Math.sin(radlat2) + Math.cos(radlat1) * Math.cos(radlat2) * Math.cos(radtheta);
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);
}
}
Upvotes: 2
Reputation: 5055
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
Upvotes: 1
Reputation: 101
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
Radian_factor <- 0.0174533
lat_1 <- (90-lat1)*Radian_factor
lat_2 <- (90-lat2)*Radian_factor
diff_long <-(lon1-lon2)*Radian_factor
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).
Upvotes: 4
Reputation: 1574
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
https://ipython.readthedocs.io/en/stable/interactive/magics.html
"""
Radius = 6.371e6
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))
s12 = Radius * c
# 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)
Upvotes: 2
Reputation: 651
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.
Upvotes: 1
Reputation: 12497
For those looking for an Excel formula based on WGS-84 & GRS-80 standards:
=ACOS(COS(RADIANS(90-Lat1))*COS(RADIANS(90-Lat2))+SIN(RADIANS(90-Lat1))*SIN(RADIANS(90-Lat2))*COS(RADIANS(Long1-Long2)))*6371
Upvotes: 3
Reputation: 964
You can calculate it by using Haversine formula which is:
a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)
c = 2 ⋅ atan2( √a, √(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
var φ1 = 28.7041.toRadians();
var φ2 = 51.5074.toRadians();
var Δφ = (51.5074-28.7041).toRadians();
var Δλ = (0.1278-77.1025).toRadians();
var a = Math.sin(Δφ/2) * Math.sin(Δφ/2) +
Math.cos(φ1) * Math.cos(φ2) *
Math.sin(Δλ/2) * Math.sin(Δλ/2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c; // metres
d = d/1000; // km
Upvotes: 1
Reputation: 5268
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));
}
Upvotes: 1
Reputation: 8979
Here is the SQL Implementation to calculate the distance in km,
SELECT UserId, ( 3959 * acos( cos( radians( your latitude here ) ) * cos( radians(latitude) ) *
cos( radians(longitude) - radians( your longitude here ) ) + sin( radians( your latitude here ) ) *
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
Upvotes: 6
Reputation: 2895
Here is a C# Implementation:
static class DistanceAlgorithm
{
const double PIx = 3.141592653589793;
const double RADIUS = 6378.16;
/// <summary>
/// Convert degrees to Radians
/// </summary>
/// <param name="x">Degrees</param>
/// <returns>The equivalent in radians</returns>
public static double Radians(double x)
{
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));
return angle * RADIUS;
}
}
Upvotes: 79
Reputation: 2139
FSharp version, using miles:
let radialDistanceHaversine location1 location2 : float =
let degreeToRadian degrees = degrees * System.Math.PI / 180.0
let earthRadius = 3959.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
* earthRadius
Upvotes: 1
Reputation: 9279
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
google.maps.geometry.spherical.computeDistanceBetween(): 900 km
geoDistance(50,5,-58,-3)
Haversine: 12030 km
Maymenn: 11135 km
Keerthana: 10310 km
google.maps.geometry.spherical.computeDistanceBetween(): 12044 km
geoDistance(.05,.005,.058,.003)
Haversine: 0.9169 km
Maymenn: 0.851723 km
Keerthana: 0.917964 km
google.maps.geometry.spherical.computeDistanceBetween(): 0.917964 km
geoDistance(.05,80,.058,80.3)
Haversine: 33.37 km
Maymenn: 33.34 km
Keerthana: 33.40767 km
google.maps.geometry.spherical.computeDistanceBetween(): 33.40770 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 R1 = radius(lat1);
var x1 = R1*cos(lat1)*cos(lng1);
var y1 = R1*cos(lat1)*sin(lng1);
var z1 = R1*sin(lat1);
var R2 = radius(lat2);
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));
}
Upvotes: 9
Reputation: 3615
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 }
radius = 6371 # Radius of the earth in kilometer
dLat = radians[to[0]-from[0]]
dLon = radians[to[1]-from[1]]
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
Upvotes: 2
Reputation: 9554
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.
For more info: wikipedia geographic distance
Upvotes: 43
Reputation: 179
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 = sin(deg2rad($lat1)) * sin(deg2rad($lat2)) + cos(deg2rad($lat1)) * cos(deg2rad($lat2)) * cos(deg2rad($theta));
$dist = acos($dist);
$dist = rad2deg($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;
}
}
Upvotes: 5
Reputation: 143
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(
sin(radians(blat-alat)/2)^2
+ (
(sin(radians(blng-alng)/2)^2) *
cos(radians(alat)) *
cos(radians(blat))
)
)
) * 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;
Upvotes: 4
Reputation: 965
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.cos(Math.toRadians(userLat)) * Math.cos(Math.toRadians(venueLat))
* Math.sin(lngDistance / 2) * Math.sin(lngDistance / 2);
double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
return (int) (Math.round(AVERAGE_RADIUS_OF_EARTH_KM * c));
}
Note that here we are rounding the answer to the nearest km.
Upvotes: 73
Related Questions
- How to find distance from the latitude and longitude of two locations?
- Calculate distance between latitude and longitude points?
- Geographic coordinate distance calculation with haversine formula gives the wrong output
- GPS - Is the haversine formula accurate for distance between two nearby gps points?
- Getting distance from longitude and latitude using Haversine's distance formula
- How to calculate the distance between two latitude and longitude
- Issues with a result from calculating latitude/longitude from haversine formula in C++
- calculate distance of two cities using Haversine formula-how to deal with minus longitudes
- Calculate distance between 2 set of lon and lat
- Distance between lat/long points using the haversine formula