Sunday 15 October 2017

inaccurate distance measurements in Python


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:


enter image description here


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

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...