一尘不染

Geod ValueError:未定义反短程线

python

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

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)

此错误消息从何而来?


阅读 160

收藏
2021-01-20

共1个答案

一尘不染

这两点仅相距几厘米。看起来pyproj/Geod不能很好地解决彼此靠近的点。这有点奇怪,因为在这样的距离下简单的​​平面几何体已经足够了。同样,该错误消息也有点可疑,因为这表明这两点是对立的,即在直径上是相反的,显然不是这种情况!OTOH,也许它所提到的对立点是在计算中以某种方式出现的中间点…不过,我还是会犹豫使用这种行为的库。

鉴于此缺陷,我怀疑pyproj还有其他缺陷。特别是,它可能会使用旧的Vincenty公式进行椭球测地线计算,众所周知,当处理近对映点时该算法不稳定,并且在大距离上并非特别精确。我建议使用CFF
Karney的现代算法。

Karney博士是一个主要贡献者上测地线的维基百科文章,特别是对一个椭球大地测量,他的geographiclib可PyPI上,这样你就可以很容易地使用安装它pip。请参阅他的SourceForge网站以获取更多信息以及其他语言的geolib绑定。

FWIW,这是一个使用geoliblib计算问题距离的简短演示。

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毫米多一点。


如果您希望看到geoliblib用于解决复杂的测地线问题,请参阅我去年在gist上使用Python
2/3源代码编写的math.stackexchange答案


希望这不再是问题,因为pyproj现在使用来自geoliblib的代码

2021-01-20