Geod ValueError : undefined inverse geodesic
Geod ValueError : undefined inverse geodesic
我想使用 pyproj
库中的 Geod
class 来计算两个经度/纬度点之间的距离。
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,即截然相反,但显然不是这样! OTOH,也许它提到的对映点是在计算中以某种方式出现的某个中间点......不过,我还是很犹豫要不要使用这样的库。
鉴于此缺陷,我怀疑 pyproj
还有其他缺陷。特别是,它可能使用旧的 Vincenty's formulae 进行椭球测地线计算,众所周知,这种方法在处理近对映点时不稳定,并且在大距离上不是特别准确。我建议使用 C. F. F. Karney 的现代算法。
Karney 博士是关于测地线的维基百科文章的主要贡献者,特别是 Geodesics on an ellipsoid, and his geographiclib is available on PyPi, so you can easily install it using pip
. See his SourceForge 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 I wrote last year, with Python 2 / 3 source code on gist。
希望这不再是问题,因为 pyproj now uses code from geographiclib。
我想使用 pyproj
库中的 Geod
class 来计算两个经度/纬度点之间的距离。
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,即截然相反,但显然不是这样! OTOH,也许它提到的对映点是在计算中以某种方式出现的某个中间点......不过,我还是很犹豫要不要使用这样的库。
鉴于此缺陷,我怀疑 pyproj
还有其他缺陷。特别是,它可能使用旧的 Vincenty's formulae 进行椭球测地线计算,众所周知,这种方法在处理近对映点时不稳定,并且在大距离上不是特别准确。我建议使用 C. F. F. Karney 的现代算法。
Karney 博士是关于测地线的维基百科文章的主要贡献者,特别是 Geodesics on an ellipsoid, and his geographiclib is available on PyPi, so you can easily install it using pip
. See his SourceForge 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 I wrote last year, with Python 2 / 3 source code on gist。
希望这不再是问题,因为 pyproj now uses code from geographiclib。