# differences in distance between shapely and geopy

I am getting totally different distances from shapely and geopy. On one hand - it makes total sense, shapely measures on a plane and geopy measures on a spheroid. But the differences are so dramatic, I am wondering if something is truly off here.

Using:

• Shapely 1.3.2
• geopy 1.9.1
• Fiona 1.l.5

my prj (for shapely/fiona) IS Web-Merc EPSG: 3857 (w/o flattening):

``````     >>print input_crs
>>{u'a': 6378137, u'lon_0': 0, u'k': 1, u'y_0': 0, u'b': 6378137, u'proj': u'merc', u'x_0': 0, u'units': u'm', u'no_defs': True}
>>print to_string(input_crs)
>>+a=6378137 +b=6378137 +k=1 +lon_0 +no_defs +proj=merc +units=m +x_0 +y_0
``````

Test Pts (*Arizona/N.Mex general area):

``````city: POINT (-12352678.6153664 3384049.040913342)
capital: POINT (-12318614.85118365 3417715.06365641)
``````

My code, using shapely for .buffer and Vincenty for the ellipsoid measure (geopy)

``````def evaluate_point(city_pt, capital_pts):

#tests whether the capitals are within a certain distance of a
#city point

# set the ellipsoid for measurement
distance.VincentyDistance.ELLIPSOID = 'WGS-84'
d = distance.distance

# create a wkt obj of city pt
city_wkt = dumps(city_pt)

dist = 300000

# create a shapely obj of all the capitals
for pt in capital_pts:
# create a shapely obj of all the captial_pts
capital_pt = shape(pt['geometry'])
# create a wkt obj of the admin1_pt
capital_to_measure = dumps(capital_pt)

## BUFFER USING SHAPELY - DIST = 300000 meters due to prj
capital_buffer = capital_pt.buffer(dist)
# city_pt is a shapely obj
if city_pt.within(capital_buffer):
distance_between_pts = d(capital_to_measure, city_wkt).meters
print"within the buffer" + " " + str(distance_between_pts)
print "city: " + city_wkt
print "capital: " + capital_to_measure

else:
print "outside the buffer"
``````

FOR this data I get this:

``````city: POINT (-12352678.6153664 3384049.040913342)
capital: POINT (-12318614.85118365 3417715.06365641)

within the buffer # shapely buffer set to (300000 meters)
14815370.7894 # Vincenty ellipsoid measured distance reported in meters
``````

## 1 Answer

1) For calculate the Vincenty distance, you need coordinates in longitude, latitude (degrees) and you choose WGS84 (look at with Geopy:distance)

2) The coordinates of your points are in meters: u'units': u'm' (according to your Proj4 string (+a=6378137 +b=6378137 +k=1 +lon_0 +no_defs +proj=merc +units=m +x_0 +y_0)

As your coordinates are in meters, you can use the Shapely distance (Euclidean Distance or Linear distance) between two points (What is the unit the shapely length attribute?)

``````for pt in capital_pts:
capital_pt = shape(pt['geometry'])
capital_to_measure = capital_pt.wkt
capital_buffer = capital_pt.buffer(dist)
if city_pt.within(capital_buffer):
distance_between_pts =  capital_to_measure.distance(city_pt) # in meters
``````

3) If you really want to use the Geopy Vincenty distance you need to change the projection of your data (with Pyproj, for example)

``````from pyproj import Proj, transform
your_proj = Proj('+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ')
wgs84 = Proj(init="epsg:4326")
city = Point(-12352678.6153664,3384049.040913342) # shapely point
long, lat =  transform(your_proj, wgs84 ,city.x, city.y)
print(Point(long,lat).wkt)
'POINT (-110.966 29.22957857026954)'
``````