我想使用Geod库中的pyproj类计算两个lon/lat点之间的距离。

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])

我得到以下错误:
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)

此错误消息从何而来?

最佳答案

这两点只有几厘米的距离。看起来pyproj/Geod不能很好地处理那些紧密相连的点。这有点奇怪,因为简单的平面几何在这样的距离上已经足够了。另外,这个错误信息有点可疑,因为它表明这两点是antipodal,即截然相反,这显然不是事实!哦,可能它提到的对极点是在计算中出现的中间点。。。不过,我在使用这样的库时还是比较犹豫的。
考虑到这个缺陷,我怀疑pyproj还有其他缺陷。特别是,它可能使用旧的Vincenty's formulae进行椭球测地线计算,这在处理近对极点时是不稳定的,在大距离上也不是特别精确。我推荐使用C.F.F.Karney的现代算法。
Karney博士是Wikipedia关于测地线的文章的主要撰稿人,特别是Geodesics on an ellipsoid,他的geographiclib可以在PyPi上找到,因此您可以使用pip轻松安装它。有关更多信息,请参见hisSourceForge site和其他语言的geographiclib绑定。
FWIW,这是一个使用geographiclib计算问题距离的简短演示。

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'])

输出
0.07345528623159624

这个数字是以米为单位的,所以这两点相距有点超过73毫米。
如果你想看到geographiclib被用来解决一个复杂的测地线问题,请看我去年写的这个math.stackexchange answer,上面有Python 2/3源代码。
希望这不再是一个问题,因为gist

08-03 20:08