Shapely Distance different from Geopy/Haversine

Question

I am running a specific analysis where we use shapely to create buffers around points (store locations) and then check if other points (user locations) are present within that buffer value. When i check the distance using shapely, it turns out to be different from the distance I get from geopy. The haversine formula agrees with Geopy and a check on google maps using the measure distance function also gives around the same distance

Here is an example:

from shapely.geometry import Point, shape
from pyproj import Proj, transform
from geopy.distance import vincenty, great_circle

pt_store=Point(transform(Proj(init='EPSG:4326'),Proj(init='EPSG:3857'),-76.799614, 39.435307))

pt_user=Point(transform(Proj(init='EPSG:4326'),Proj(init='EPSG:3857'),-76.79989,39.43604))

vincenty((39.435307,-76.799614),(39.43604,-76.79989)).meters
great_circle((39.435307,-76.799614),(39.43604,-76.79989)).meters
pt_store.distance(pt_user)

Vincenty: 84.77847691521336
Great_circle: 84.90640111682812
Shapely: 110.02637304449682
Haversine formula (http://www.movable-type.co.uk/scripts/latlong.html): 84.88

Which one is right? Shapely or others? Also, is such a big difference(~22%) expected? Or am I missing something?

Answer 1

Because the principles and the algorithms are different (look at Geographical distance)

1. Shapely use the euclidean distance in a cartesian plane and the shortest distance between two points in a plane is a straight line which contains the two points.

 import numpy as np
 print np.linalg.norm(np.array(pt_user) - np.array(pt_store))
 110.02637304449682 # meters
 from scipy.spatial import distance
 print distance.euclidean(pt_user, pt_store)
 110.02637304449682 # meters

2. Vincenty, Great Circle and Haversine use either the geodesic distance (on an ellipsoid, Vincenty) or the great-circle distance (the shortest distance along the surface of a sphere) between two points. The shortest distance on the surface of a sphere is along the great-circle which contains the two points.

   Therefore it is normal that the Shapely, Numpy and Scipy euclidean distances differ from the Vincenty, Great Circle and Haversine distances and the differences between the Vincenty, Great Circles and Haversine distances are linked to the choice of an ellipsoid, and many other things.

   You can also change the ellipsoid

 print vincenty((39.435307,-76.799614),(39.43604,-76.79989),ellipsoid='WGS-84')
 0.0847784769149 km
 print vincenty((39.435307,-76.799614),(39.43604,-76.79989),ellipsoid='GRS-80')
 0.0847784769128 km

Or use other libraries as geodistance

 print geodistance.distanceVincenty(39.435307,-76.799614,39.43604,-76.79989, ellipsoid='WGS-84')
 (0.08477847691523362, -16.276730447136675) # distance, azimuth
 print geodistance.distanceHaversine(39.435307,-76.799614,39.43604,-76.79989)
 (0.08488248586585143, -16.214988211007256)

You can see that all the differences are centimetric. With metric precision, all the values = 85 meters.

3. Which one is right? All, because it depends on the context: if you work with projected data (cartesian plane), you use the Euclidean distance (Shapely, Numpy ,Scipy and many others), if not, one of the others.

   They are also many other distances (Scipy Spatial distances)

New

In support of the answer of Mintx

pt_store=Point(transform(Proj(init='EPSG:4326'),Proj(init='EPSG:31370'),-76.799614, 39.435307))
pt_user=Point(transform(Proj(init='EPSG:4326'),Proj(init='EPSG:31370'),-76.79989,39.43604))
pt_store.distance(pt_user)
86.26511001003892

Answer 2

And here's another distance calculation from GeographicLib:

from geographiclib.geodesic import Geodesic
g = Geodesic.WGS84.Inverse(39.435307, -76.799614, 39.43604, -76.79989)
print(g['s12'])  # 84.7784769689

I would consider this the right one, within 15 nanometres.

Answer 3

When implementing haversine I always trip on this, changing latitude and longitude in the implementation. Here is alll the three

from geopy import distance
cord1 = (-81.81269,39.33015)
cord2 = (-80.4353 ,37.6465)
distance_m = distance.geodesic(cord1,cord2).m
print(f"The Total distance from from geopy.geodesic is {distance_m/1000:.2f} km")

distance_m = distance.great_circle(cord1,cord2).m
print(f"The Total distance from from geopy.great_circle is {distance_m/1000:.2f} km")

def haversine(lon1: float, lat1: float, lon2: float, lat2: float) -> float:
    RADIUS_OF_EARTH_IN_M = 6371.01 *1000
    """
    Calculate the distance between two points on the Earth's surface using the Haversine formula.

    Args:
        lon1 (float): Longitude of the first point in degrees.
        lat1 (float): Latitude of the first point in degrees.
        lon2 (float): Longitude of the second point in degrees.
        lat2 (float): Latitude of the second point in degrees.

    Returns:
        float: The distance between the two points in meters.

    Notes:
        The Haversine formula calculates the great-circle distance between two points on a sphere, such as the Earth.
        The calculated distance is an approximation and assumes a perfect sphere.

    """
    dLat = np.radians(lat2 - lat1)
    dLon = np.radians(lon2 - lon1)
    lat1 = np.radians(lat1)
    lat2 = np.radians(lat2)
 
    a = np.sin(dLat/2)**2 + np.cos(lat1)*np.cos(lat2)*np.sin(dLon/2)**2
    c = 2*np.arcsin(np.sqrt(a))
 
    return RADIUS_OF_EARTH_IN_M * c


distance_m  =haversine(cord1[1],cord1[0],cord2[1],cord2[0])
print(f"The Total distance from  haversine is {distance_m/1000:.2f} km")

Output

The Total distance from from geopy.geodesic is 156.51 km
The Total distance from from geopy.great_circle is 155.84 km
The Total distance from  haversine is 155.84 km