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
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)
No comments:
Post a Comment