ST_Distance_Sphere() in Python?

Question

I am working on a Python project where I have two lat/long pairs and I want to calculate the distance between them. In other projects I have calculated distance in Postgres using ST_Distance_Sphere(a.loc_point, b.loc_point), but I would like to avoid having to load all of my data into Postgres just so that I can calculate distance differences. I have searched, but have not been able to find what I would like, which is a purely Python implementation of this so that I don't have to load my data into Postgres.

I know there are other distance calculations that treat the earth as a perfect sphere, but those aren't good enough due to poor accuracy, which is why I would like to use the PostGIS ST_Distance_Sphere() function (or an equivalent).

Here are a couple of sample Lat/Longs that I would like to calculate the distance of:

Lat, Long 1: (49.8755, 6.07594)
Lat, Long 2: (49.87257, 6.0784)

I can't imagine I am the first person to ask this, but does anyone know of a way to use ST_Distance_Sphere() for lat/long distance calculations purely from within a Python script?

ewcz · Answer 1 · 2017-02-16T15:53:28.773

3

I would recommend the geopy package - see section Measuring Distance in the documentation...

For your particular case:

from geopy.distance import great_circle

p1 = (49.8755, 6.07594)
p2 = (49.87257, 6.0784)

print(great_circle(p1, p2).kilometers)

edited Feb 16 '17 at 15:53

answered Feb 16 '17 at 15:42

ewcz

12,819
1
25
47

Sishaar Rao · Accepted Answer · 2017-02-16T15:20:24.203

1

This is a rudimentary function used to calculate distance between two coordinates on a perfect sphere with Radius = Radius of Earth

from math import pi , acos , sin , cos
def calcd(y1,x1, y2,x2):
   #
   y1  = float(y1)
   x1  = float(x1)
   y2  = float(y2)
   x2  = float(x2)
   #
   R   = 3958.76 # miles
   #
   y1 *= pi/180.0
   x1 *= pi/180.0
   y2 *= pi/180.0
   x2 *= pi/180.0
   #
   # approximate great circle distance with law of cosines
   #
   x = sin(y1)*sin(y2) + cos(y1)*cos(y2)*cos(x2-x1)
   if x > 1:
       x = 1
   return acos( x ) * R

Hope this helps!

edited Feb 16 '17 at 15:20

answered Feb 16 '17 at 15:08

Sishaar Rao

340
1
8

How can you get answer in meters with this solution? – Chris Nielsen Feb 16 '17 at 15:57
In the line R = 3958.76 # miles simply set R to the radius of the earth in meters (6.371 million meters) and that should return answers in terms of meters – Sishaar Rao Feb 16 '17 at 15:59

score 1 · Answer 3 · edited May 23 '17 at 11:53

See this How can I quickly estimate the distance between two (latitude, longitude) points?

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)) 
    km = 6367 * c
    return km

By Aaron D

You can modify it to return miles by adding miles = km * 0.621371

Chris Nielsen · Answer 4 · 2017-02-16T20:05:23.000

I have since found another way in addition to the answers provided here. Using the python haversine module.

from haversine import haversine as h

# Return results in meters (*1000)
print '{0:30}{1:12}'.format("haversine module:", h(a, b)*1000)

I tested all three answers plus haversine module against what I got using ST_Distance_Sphere(a, b) in Postgres. All answers were excellent (thank you), but the all math answer (calcd) from Sishaar Rao was the closest. Here are the results:

# Short Distance Test
ST_Distance_Sphere(a, b):     370.43790478    
vincenty:                     370.778186438
great_circle:                 370.541763803
calcd:                        370.437386736
haversine function:           370.20481753
haversine module:             370.437394767

#Long Distance test:
ST_Distance_Sphere(a, b):     1011734.50495159
vincenty:                     1013450.40832
great_circle:                 1012018.16318
calcd:                        1011733.11203
haversine function:           1011097.90053
haversine module:             1011733.11203

ST_Distance_Sphere() in Python?

4 Answers4