66

I want to be able to get a estimate of the distance between two (latitude, longitude) points. I want to undershoot, as this will be for A* graph search and I want it to be fast. The points will be at most 800 km apart.

fread2281
  • 1,136
  • 1
  • 11
  • 31
  • 1
    Should we infer these points lie on a *sphere*? – phs Apr 01 '13 at 02:45
  • 1
    See http://stackoverflow.com/questions/27928/how-do-i-calculate-distance-between-two-latitude-longitude-points or http://stackoverflow.com/questions/4913349/haversine-formula-in-python-bearing-and-distance-between-two-gps-points (python) – Aaron D Apr 01 '13 at 02:45
  • Yes, on earth, but speed. AFAIK complex math is not fast enough. – fread2281 Apr 01 '13 at 02:48
  • 14
    I suggest you measure first before concluding it's not fast enough. – phs Apr 01 '13 at 02:50
  • Ok, but also my latitude, longitude is not 100% accurate, so I need to underestimate a bit. – fread2281 Apr 01 '13 at 02:56
  • 1
    Sometimes it's possible to know enough about an implementation and algorithm to know performance won't be good enough even prior to benchmarking. For instance, one case where the haversine distance method isn't appropriate is when attempting to match large datasets on proximity, as the haversine algorithm doesn't allow any predicate pushdowns or partition matching in most querying engines. We found that leveraging approximate distances with pushdowns to produce a cartesian clustering base took ~1/50th the time on a 250k record dataset. The accepted answer would take over a week to run here. – bsplosion Mar 13 '19 at 13:18
  • As I mention [here](https://stackoverflow.com/a/71412448/7109869), one can use the Haversine Formula. However, keep in mind that one is approximating the Earth as a sphere, and that has errors ([see this answer](https://gis.stackexchange.com/a/25580/130105)) – Gonçalo Peres Mar 09 '22 at 16:13

8 Answers8

124

The answers to Haversine Formula in Python (Bearing and Distance between two GPS points) provide Python implementations that answer your question.

Using the implementation below I performed 100,000 iterations in less than 1 second on an older laptop. I think for your purposes this should be sufficient. However, you should profile anything before you optimize for performance. 

from math import radians, cos, sin, asin, sqrt
def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance between two points 
    on the earth (specified in decimal degrees)
    """
    # convert decimal degrees to radians 
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a)) 
    # Radius of earth in kilometers is 6371
    km = 6371* c
    return km

To underestimate haversine(lat1, long1, lat2, long2) * 0.90 or whatever factor you want. I don't see how introducing error to your underestimation is useful.

RoyaumeIX
  • 1,947
  • 4
  • 13
  • 37
Aaron D
  • 5,817
  • 1
  • 36
  • 51
  • So, this isn't a "quick" method, it's the "standard" method, correct? (and slower than my equally accurate method, when comparing in Excel anyhow) – ashleedawg Sep 20 '18 at 21:53
  • 2
    @ashleedawg - The Haversine formula is a particular method for estimating distance between two points assuming Earth is a sphere. I invite you to look at the wikipedia page on Haversine formula if you want to explore other methods. Haversine is common but there is not a "standard" formula. If your "equally accurate method" is good, why not post it as an answer here? – Aaron D Oct 15 '18 at 22:40
45

Since the distance is relatively small, you can use the equirectangular distance approximation. This approximation is faster than using the Haversine formula. So, to get the distance from your reference point (lat1, lon1) to the point you're testing (lat2, lon2) use the formula below:

from math import sqrt, cos, radians

R = 6371  # radius of the earth in km
x = (radians(lon2) - radians(lon1)) * cos(0.5 * (radians(lat2) + radians(lat1)))
y = radians(lat2) - radians(lat1)
d = R * sqrt(x*x + y*y)

Since R is in km, the distance d will be in km.

Reference: http://www.movable-type.co.uk/scripts/latlong.html

wim
  • 338,267
  • 99
  • 616
  • 750
TreyA
  • 3,339
  • 2
  • 20
  • 26
  • 1
    Note that the equirectangular approximation will overshoot, not undershoot. You can derive it from the Haversine formula using small angle approximations ``sin^2(dlon) ~ dlon^2``, ``sin^2(dlat) ~ dlat^2`` and ``cos(dlat) ~ 1`` where ``dlon=lon2-lon1`` and ``dlat=lat2-lat1``. All approximations are >= the exact version, thus the approximated distance will be larger than the exact distance. – sfinkens Aug 01 '18 at 07:30
  • I could argue that Haversine (if that is what you mean by 'exact') is an approximation as well since the Earth isn't a true sphere. And what order of magnitude of overshoot are we talking? A couple of inches overshoot? feet? – TreyA Aug 01 '18 at 17:05
  • Of course, that is correct. The OP asked for undershooting. But for some applications overshooting is even a desired property of the approximation, e.g. for pre-selecting neighbouring points within a given distance. So the point is that it overshoots, not how much. – sfinkens Aug 02 '18 at 07:55
  • Yes, but if it overshoots by only inches, technically you are correct but practically who cares? The OP is looking for speed as well. And you are free to adjust the radius of the Earth to your under/over shooting liking. – TreyA Aug 02 '18 at 15:06
8

One idea for speed is to transform the long/lat coordinated into 3D (x,y,z) coordinates. After preprocessing the points, use the Euclidean distance between the points as a quickly computed undershoot of the actual distance.

Raymond Hettinger
  • 216,523
  • 63
  • 388
  • 485
7

If the distance between points is relatively small (meters to few km range) then one of the fast approaches could be

from math import cos, sqrt
def qick_distance(Lat1, Long1, Lat2, Long2):
    x = Lat2 - Lat1
    y = (Long2 - Long1) * cos((Lat2 + Lat1)*0.00872664626)  
    return 111.319 * sqrt(x*x + y*y)

Lat, Long are in radians, distance in km.

Deviation from Haversine distance is in the order of 1%, while the speed gain is more than ~10x.

0.00872664626 = 0.5 * pi/180,

111.319 - is the distance that corresponds to 1degree at Equator, you could replace it with your median value like here https://www.cartographyunchained.com/cgsta1/ or replace it with a simple lookup table.

Danylo Zherebetskyy
  • 1,429
  • 13
  • 9
  • Can someone please clarify how to get this magic 111.138 number? Distance of 1degree at Equator. How would I get the appropriate value for somewhere else? Found this, maybe helpful... https://blog.mapbox.com/fast-geodesic-approximations-with-cheap-ruler-106f229ad016 – Liam Mitchell Jul 25 '19 at 08:48
  • 1
    @LiamMitchell thank you for noticing the typo, I've corrected the Equatorial distance: Equatorial radius (km)= 6378.137 accordingly to https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html, hence 2*pi*R / 360 = 111.319491 km – Danylo Zherebetskyy Sep 20 '19 at 18:53
3

For maximal speed, you could create something like a rainbow table for coordinate distances. It sounds like you already know the area that you are working with, so it seems like pre-computing them might be feasible. Then, you could load the nearest combination and just use that.

For example, in the continental United States, the longitude is a 55 degree span and latitude is 20, which would be 1100 whole number points. The distance between all the possible combinations is a handshake problem which is answered by (n-1)(n)/2 or about 600k combinations. That seems pretty feasible to store and retrieve. If you provide more information about your requirements, I could be more specific.

Aaron D
  • 5,817
  • 1
  • 36
  • 51
1

You can use cdist from scipy spacial distance class:

For example:

from scipy.spatial.distance import cdist 
df1_latlon = df1[['lat','lon']]
df2_latlon = df2[['lat', 'lon']]
distanceCalc = cdist(df1_latlon, df2_latlon, metric=haversine)
Grayrigel
  • 3,474
  • 5
  • 14
  • 32
Ishaan Sharma
  • 55
  • 1
  • 8
  • 2
    As far as I see in the cdist doc (docs.scipy.org/doc/scipy/reference/generated/…), there is no haversine metric defined. Running it also confirms that by "Unknown Distance Metric: haversine" Error. – pegah Aug 16 '21 at 15:47
  • Correct. Just cross-checked . Only following are available : braycurtis,canberra,chebyshev,cityblock,correlation,cosine,euclidean,jensenshannon,mahalanobis,minkowski,seuclidean,sqeuclidean. https://github.com/scipy/scipy/blob/master/scipy/spatial/distance.py – Ishaan Sharma Dec 16 '21 at 04:48
  • instead you can use this for finding the haversine distance : https://github.com/mapado/haversine#calculate-the-distance-between-lyon-and-paris – Ishaan Sharma Dec 16 '21 at 04:53
0

To calculate a haversine distance between 2 points u can simply use mpu.haversine_distance() library, like this: 

>>> import mpu
>>> munich = (48.1372, 11.5756)
>>> berlin = (52.5186, 13.4083)
>>> round(mpu.haversine_distance(munich, berlin), 1)
>>> 504.2
Chiefir
  • 2,561
  • 1
  • 27
  • 46
-1

Please use the following code.

def distance(lat1, lng1, lat2, lng2):
    #return distance as meter if you want km distance, remove "* 1000"
    radius = 6371 * 1000 

    dLat = (lat2-lat1) * math.pi / 180
    dLng = (lng2-lng1) * math.pi / 180

    lat1 = lat1 * math.pi / 180
    lat2 = lat2 * math.pi / 180

    val = sin(dLat/2) * sin(dLat/2) + sin(dLng/2) * sin(dLng/2) * cos(lat1) * cos(lat2)    
    ang = 2 * atan2(sqrt(val), sqrt(1-val))
    return radius * ang
uher
  • 19
  • 2
  • In my case, the other codes is not working well for me. So, I just transtlate the funtion in this answer http://stackoverflow.com/questions/6981916/how-to-calculate-distance-between-two-locations-using-their-longitude-and-latitu – uher Aug 02 '16 at 07:28