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.

  • The UTM zone is easy to "work out" based on the longitudes. Pick a typical longitude lambda, -180 <= lambda < 180, and use it to compute the zone number as Int((180+lambda)/6)+1. Use the sign of the latitude to decide between north and south. You don't need to use the special polar zones in high latitudes; in fact, really close to a pole you can use almost any UTM zone. – whuber Mar 15 '11 at 21:35
up vote 17 down vote accepted

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.

  • That sounds pretty perfect, I guess I need to make my own init string for proj4 then, rather than being able to use any of the existing EPSG strings? – Karl P Mar 15 '11 at 18:01
  • 1
    +1 adapt projection to the points. I prefer transverse plate carrée over transverse Mercator, but over small enough areas ("large scale"), almost any projection "centered" near the region of interest will give good accuracy. – David Cary Mar 15 '11 at 18:27
  • @David Interesting idea. On the sphere the transverse Plate Carree (Cassini) will be close to the approximate formula I gave at gis.stackexchange.com/posts/2964/edit (which might be an acceptable solution here). Formulas for TM and TPC are similar on the sphere; on an ellipsoid, TPC is slightly simpler. TM is probably supported by more software. – whuber Mar 15 '11 at 19:04
  • @Karl You could use any of the TM zones if you wish. Just shift all longitudes so that a central point in your region of interest coincides with the central meridian of the chosen zone. Multiply all distances by 1/0.9996 (and multiply all areas by the square of this factor), don't change any angles or bearings, and--if your calculations produce new points--just shift their longitudes back to the original coordinate system. – whuber Mar 15 '11 at 19:07

I haven't tried this but from the documentation it looks like you can use http://search.cpan.org/~grahamc/Geo-Coordinates-UTM-0.08/UTM.pm#latlon_to_utm to get from a lat/lon pair (plus ellipsoid) to UTM Zone and coordinates list. Then you can carry on with your calculation as before.

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.