python - Geod ValueError : undefined inverse geodesic

Question

I want to compute the distance between two lon / lat points by using Geod class from pyproj library.

from pyproj import Geod

g = Geod(ellps='WGS84')
lonlat1 = 10.65583081724002, -7.313341167341917
lonlat2 = 10.655830383300781, -7.313340663909912

_, _, dist = g.inv(lonlat1[0], lonlat1[1], lonlat2[0], lonlat2[1])

I get the following error :

ValueError                                Traceback (most recent call last)
<ipython-input-5-8ba490aa5fcc> in <module>()
----> 1 _, _, dist = g.inv(lonlat1[0], lonlat1[1], lonlat2[0], lonlat2[1])

/usr/lib/python2.7/dist-packages/pyproj/__init__.pyc in inv(self, lons1, lats1, lons2, lats2, radians)
    558         ind, disfloat, dislist, distuple = _copytobuffer(lats2)
    559         # call geod_inv function. inputs modified in place.
--> 560         _Geod._inv(self, inx, iny, inz, ind, radians=radians)
    561         # if inputs were lists, tuples or floats, convert back.
    562         outx = _convertback(xisfloat,xislist,xistuple,inx)

_geod.pyx in _geod.Geod._inv (_geod.c:1883)()

ValueError: undefined inverse geodesic (may be an antipodal point)

Where does this error message come from ?

Answer 1

Those two points are only a few centimetres apart. It looks like pyproj / Geod doesn't cope well with points which are that close together. That's a bit strange, since simple plane geometry is more than adequate at such distances. Also, that error message is a bit suspicious, since it's suggesting that the two points are antipodal, i.e., diametrically opposite, which is clearly not the case! OTOH, maybe the antipodal point it mentions is some intermediate point that arises somehow in the calculation... Still, I'd be rather hesitant in using a library that behaves like this.

Given this defect, I suspect that pyproj has other flaws. In particular, it probably uses the old Vincenty's formulae for its ellipsoid geodesic calculations, which is known to be unstable when dealing with near-antipodal points, and not particularly accurate over large distances. I recommend using the modern algorithms of C. F. F. Karney.

Dr Karney is a major contributor to the Wikipedia articles on geodesics, in particular Geodesics on an ellipsoid, and his geographiclib is available on PyPi, so you can easily install it using pip. See his SourceForge site for further information, and geographiclib binding in other languages.

FWIW, here's a short demo of using geographiclib to compute the distance in your question.

from geographiclib.geodesic import Geodesic

Geo = Geodesic.WGS84

lat1, lon1 = -7.313341167341917, 10.65583081724002
lat2, lon2 = -7.313340663909912, 10.655830383300781

d = Geo.Inverse(lat1, lon1,  lat2, lon2)
print(d['s12'])

output

0.07345528623159624

That figure is in metres, so those two points are a little over 73mm apart.

---

If you'd like to see geographiclib being used to solve a complex geodesic problem, please see this math.stackexchange answer I wrote last year, with Python 2 / 3 source code on gist.

---

Hopefully, this is no longer an issue, since pyproj now uses code from geographiclib.