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[0], point1_geo[1])
point2 = proj(point2_geo[0], point2_geo[1])
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.
Best Answer
Yes, you will get these kinds of errors with a global Mercator projection: it is accurate at the equator and the distortion increases exponentially with latitude away from the equator. The distance distortion is exactly 2 (100%) at 60 degrees latitude. At your test latitudes (64.14 degrees) I compute a distortion of 2.294, exactly agreeing with the ratio 904/394 = 2.294. (Earlier I computed 2.301 but that was based on a sphere, not the WGS84 ellipsoid. The difference (of 0.3%) gives us a sense of the accuracy you might gain from using an ellipsoid-based projection versus the sphere-based Haversine formula.)
There is no such thing as a global projection that yields highly accurate distances everywhere. That's one reason the UTM zone system is used!
One solution is to use spherical geometry for all your calculations, but you have rejected that (which is reasonable if you're going to be doing complex operations, but the decision might be worth revisiting).
Another solution is to adapt the projection to the points being compared. For example, you could safely use a transverse Mercator (as in the UTM system) with a meridian lying near the center of the region of interest. Moving the meridian is a simple thing to do: just subtract the meridian's longitude from all the longitudes and use a single TM projection centered at the Prime Meridian (with a scale factor of 1, rather than the 0.9996 of the UTM system). For your work this will tend to be more accurate than using UTM itself. It will give correct angles (TM is conformal) and will be remarkably accurate for points separated by only a few tens of kilometers: expect better than six-digit accuracy. In fact, I would be inclined to attribute any small differences between these adapted-TM distances and the Haversine distances to the difference between the ellipsoid (used for the TM projection) and the sphere (used by Haversine), rather than to distortion in the projection.