diferencias de distancia entre shapely y geopy

Question

Estoy obteniendo distancias totalmente diferentes de shapely y geopy. Por un lado tiene todo el sentido, shapely mide en un plano y geopy mide en un esferoide . Pero las diferencias son tan dramáticas que me pregunto si algo está realmente mal aquí.

Utilizando:

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

mi prj (para shapely/fiona) ES Web-Merc EPSG: 3857 (sin aplanar):

     >>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

Pts de prueba (*área general de Arizona/N.Mex):

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

Mi código, utilizando shapely para .buffer y Vincenty para la medida del elipsoide (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"

Por estos datos obtengo esto:

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

Answer 1

1) Para calcular la distancia de Vincenty, necesitas coordenadas en longitud, latitud (grados) y eliges WGS84 (mira con Geopy:distancia )

2) Las coordenadas de sus puntos están en metros: u'unidades': u'm' (según su cadena Proj4 (+a=6378137 +b=6378137 +k=1 +lon_0 +no_defs +proj=merc +units=m +x_0 +y_0)

Como sus coordenadas están en metros, puede utilizar la distancia Shapely (distancia euclidiana o distancia lineal) entre dos puntos ( ¿Cuál es la unidad del atributo de longitud de forma? )

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) Si realmente quieres utilizar la distancia Geopy Vincenty tienes que cambiar la proyección de tus datos (con Pyproj , por ejemplo)

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)'