问题描述
我想通过使用pyproj库中的Geod类来计算两个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不能很好地应对彼此靠近的点.这有点奇怪，因为在这样的距离下简单的平面几何体已经足够了.另外，该错误消息也有点可疑，因为它表明这两点是对映，即，完全相反，显然不是这种情况！ OTOH，也许它提到的对立点是在计算中以某种方式出现的中间点...不过，我还是会犹豫使用这样的库.
鉴于此缺陷，我怀疑pyproj还有其他缺陷.特别是，它可能会使用旧的 Vincenty公式进行椭球测地线计算，即已知在处理近对映点时不稳定，并且在大距离上不是特别准确.我建议使用C. F. F. Karney的现代算法.
Karney博士是Wikipedia关于测地线的文章的主要撰稿人，特别是椭球上的测地学，而他的 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用于解决复杂的测地线问题，请参见此 math.stackexchange答案去年我在要点上用Python 2/3源代码编写了... >
希望这不再是问题，因为 pyproj现在使用来自geoliblib的代码.
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 ?
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.
这篇关于Geod ValueError:未定义反短程线的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！