inaccurate distance measurements in Python

Question

I've tried several ways of measuring the distance between two points in my Django application and comparing the results to a reliable measurement. My numbers are way off. For example, I'm assuming the real distance between Dublin and Liverpool is ~217 km as reported by Google Maps:

Using geopy (should be extremely accurate):

>>> from geopy.distance import distance
>>> dublin = (-6.270447, 53.339791)
>>> liverpool = (-2.991028, 53.402061)
>>> distance(dublin, liverpool).km
362.70989455939394

Using Django's GEOS API (less accurate linear calculation):

>>> from django.contrib.gis.geos import Point
>>> dublin = Point(-6.270447, 53.339791, srid=3857)
>>> liverpool = Point(-2.991028, 53.402061, srid=3857)
>>> dublin.distance(liverpool)*100
328.00101418228877

EDIT: Using a better projection for this area (UTM 30N) yields almost the same result:

>>> dublin.transform(32630)
>>> liverpool.transform(32630)
>>> dublin.distance(liverpool)*100
328.32200116442266

In both cases I'm off by over 100km! Measuring small distances (<1km) is just as inaccurate. What am I doing wrong here!?

Answer 1

If you inverse the coordinates, it does not work (geopy uses (latitude,longitude) in the WGS84 crs)

 dublin = (53.33306,-6.24889)
 liverpool  = ( 53.41058,-2.97794)
 print distance(dublin, liverpool).km
 217.863019038
 print(vincenty(dublin, liverpool).kilometers)
 217.863019038
 print(great_circle(dublin, liverpool).kilometers)
 217.211596704

GEOS (shapely, django) uses a Cartesian plane and the Euclidean distance. With pyproj (django uses (longitude,latitude))

from django.contrib.gis.geos import Point
dublin = Point(-6.270447, 53.339791, srid=4326) # in degrees
liverpool = Point(-2.991028, 53.402061, srid=4326) # in degrees
dublin.distance(liverpool)*100
328.00101418228877 # units ?
import pyproj
# conversion from WGS84 to epsg:3857
p1 = pyproj.Proj(proj='latlong',datum='WGS84')
p2 = pyproj.Proj(init='epsg:3857')
a = pyproj.transform(p1,p2,-6.270447, 53.339791)
b = pyproj.transform(p1,p2,-2.991028, 53.402061)
dublin = Point(a) # in meters
liverpool = Point(b) # in meters
dublin.distance(liverpool)/1000 # Euclidean
365.2480859440489 #in km

But, as Vince says, the Mercator projection should not ever be used for distance measurement.

With the EPSG:32630 (UTM zone 30N):

p3 = pyproj.Proj(init='epsg:32630')
a = pyproj.transform(p1,p3,-6.270447, 53.339791)
b = pyproj.transform(p1,p3,-2.991028, 53.402061)
dublin = Point(a)
liverpool = Point(b)
dublin.distance(liverpool)/1000
218.32514783088294 #in km

And all the results (geopy and django ) are comparable with the Google distance or the Distance from Liverpool to Dublin (218 km)