[GIS] differences in distance between shapely and geopy

fionageopypyqgisqgisshapely

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

Best 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)'
Related Question