I've got points in WGS84 lat/long and I'd like to measure "small" (less than say 5km) distances between them.
I can use the haversine formula from http://www.movable-type.co.uk/scripts/latlong.html and it works very well.
I'd like to use Python Shapely libraries though, so that I can do more operations than just distance, and because at the scale I'm working with, a flat earth is a good enough approximation. To reliably project the geographic coords to a cartesian coord, I'm using Python's
proj4, but seem to get bigger errors than I'd like.
If I use the local UTM zone, I get differences between haversine of a couple of meters, which is fine. But I don't want to have to work out the UTM zone (the points could be worldwide), so I tried with "spherical Mercator" but now the differences between haversine and projected distances are well over 100%. Is this really right for spherical Mercator? All I really want is a workable Cartesian projection for two points within 5km of each other anywhere in the world.
from shapely.geometry import Point from pyproj import Proj proj = Proj(proj='utm',zone=27,ellps='WGS84') #proj = Proj(init="epsg:3785") # spherical mercator, should work anywhere... point1_geo = (-21.9309694, 64.1455718) point2_geo = (-21.9372481, 64.1478206) point1 = proj(point1_geo, point1_geo) point2 = proj(point2_geo, point2_geo) point1_cart = Point(point1) point2_cart = Point(point2) print "p1-p2 (haversine)", hdistance(point1_geo, point2_geo) print "p1-p2 (cartesian)", point1_cart.distance(point2_cart)
At this point, the haversine distance between them is 394m, and using utm zone 27, 395m. But if I use spherical Mercator, the Cartesian distance is 904m, which is way off.