如何一次计算沿路径（lat/lon 点）的测地距离？
How to calculate geodesic distance along a path (lat/lon points) at once?
在 Win10 x64 上使用 Python 3.7 和 Jupyter Notebook。
我有一个表示路径的经纬度元组列表，我想计算沿该路径的总长度（以米为单位）。我想避免计算每个段距离然后将它们加在一起，因为我必须对 700 万条路径执行此操作。所以时间效率是关键。在计算单个距离后添加所有段每条路径需要 7 毫秒。我想让它至少快 1 毫秒。
编辑： 我需要使用 WGS84 椭球来计算距离，因此球面 (Haversine) 是不够的。我想我可以以 1m 的精度工作。点沿路径随机分布。有的可以20公里，有的不到1公里。
这是带点的路径（十进制lat/lon）：
[(49.009722, 2.547778), (49.015556, 2.573611), (49.021389, 2.599167), (49.039167, 2.676389), (49.048056, 2.715), (49.044444, 2.835), (49.041667, 2.928333), (49.042778, 2.942222), (49.051667, 3.066667), (49.061389, 3.205), (49.072222, 3.357222), (49.085, 3.536944), (49.086111, 3.550833), (49.097778, 3.729444), (49.113056, 3.963056), (49.130833, 4.238056), (49.138889, 4.361667), (49.1925, 4.564444), (49.306667, 4.995556), (49.333611, 5.096944), (49.395, 5.329167), (49.490556, 5.690833), (49.514444, 5.781111), (49.53, 5.845833), (49.599444, 6.127778), (49.637222, 6.281667), (49.673333, 6.440278), (50.0475, 8.078333), (50.053611, 8.637222), (50.056667, 8.800278), (50.063056, 9.19), (50.066944, 9.486389), (50.07, 9.783056), (50.072778, 10.098611), (50.073333, 10.242778), (50.075278, 10.728889), (50.046667, 10.863333), (50.0325, 10.930278), (49.981111, 11.172222), (49.969722, 11.225833), (49.961111, 11.491389), (49.959444, 11.547222), (49.957222, 11.617222), (49.946111, 11.9325), (49.937222, 12.343889), (49.933333, 12.47), (49.9325, 12.498056), (49.928611, 12.624167), (49.924167, 12.764444), (49.919444, 12.918611), (49.910833, 13.199167), (49.909444, 13.241111), (49.907778, 13.283611), (49.900556, 13.481944), (50.077222, 13.840278), (50.124167, 13.995556), (50.182778, 14.189722), (50.220278, 14.315), (50.268889, 14.478056), (50.211389, 14.403611), (50.166389, 14.345), (50.133611, 14.3025), (50.100833, 14.26)]
我发现了cartopy - a package providing geodesic calculations that uses shapely and proj。然而，关于 cartpopy 的文档没有给出相关示例，我被困在了这一点上。基本上 geometry_length 一次性给出了一个有形物体的长度，所以我用下面的方式来做：
#defining the geoid on which to make calculations
myGeod = geodesic.Geodesic(6378137.0,1 / 298.257223563)
#making my list of latlon (in decimal degrees) into a shapely
shapelyObject = LineString(list(latlon_dd))
#applying the method on the shapelyObject given the defined ellipsoid
myGeod.geometry_length(shapelyObject)
我想计算长度（以米为单位），应该在 917,315.3 米左右。相反，我得到这个 ValueError:
ValueError Traceback (most recent call last)
<ipython-input-243-7c75042775e3> in <module>
6
7 #applying the method on the shapelyObject given the defined ellipsoid
----> 8 myGeod.geometry_length(shapelyObject)
lib\cartopy\geodesic\_geodesic.pyx in cartopy.geodesic._geodesic.Geodesic.geometry_length()
lib\cartopy\geodesic\_geodesic.pyx in cartopy.geodesic._geodesic.Geodesic.geometry_length()
lib\cartopy\geodesic\_geodesic.pyx in cartopy.geodesic._geodesic.Geodesic.inverse()
ValueError: Expecting input points to be (N, 2), got (1, 63)
提前致谢！
将根据找到的回复回答我的问题 here。显然这是一个错误，目前的解决方法是使用：
myGeod.geometry_length(np.array(shapelyObject.coords))
而不是
myGeod.geometry_length(shapelyObject)
将在最终解决方案可用时更新。
您可能不想包含另一个库，但是 LatLon 有一个内置方法来计算两个 lat/lon 点之间的 WGS84 距离（以千米为单位）。
示例取自链接页面：
>> palmyra = LatLon(Latitude(5.8833), Longitude(-162.0833)) # Location of Palmyra Atoll
>> honolulu = LatLon(Latitude(21.3), Longitude(-157.8167)) # Location of Honolulu, HI
>> distance = palmyra.distance(honolulu) # WGS84 distance in km
>> print distance
1766.69130376
>> print palmyra.distance(honolulu, ellipse = 'sphere') # FAI distance in km
1774.77188181
编辑：发帖后才注意到你想计算一条路径的距离，其中有很多点不在两个点之间（一次），抱歉...
您仍然可以使用 pyproj 并一次解决整个列表
geod = pyproj.Geod(ellps='WGS84')
_, _, distances_in_meters = geod.inv(
lons1_float_or_list_or_numpy_array,
lats1_float_or_list_or_numpy_array,
lons2_float_or_list_or_numpy_array,
lats2_float_or_list_or_numpy_array)
要从您提供的内容中获取此格式，我们可以使用 4 个理解和 ittertools.combinations。
your_list = [(49.009722, 2.547778), (49.015556, 2.573611), ...]
edges = ittertools.combinations(your_list, 2)
lons1 = [edge[0][1], for edge in edges]
lats1 = [edge[0][0], for edge in edges]
lons2 = [edge[1][1], for edge in edges]
lats2 = [edge[1][0], for edge in edges]
在 Win10 x64 上使用 Python 3.7 和 Jupyter Notebook。
我有一个表示路径的经纬度元组列表，我想计算沿该路径的总长度（以米为单位）。我想避免计算每个段距离然后将它们加在一起，因为我必须对 700 万条路径执行此操作。所以时间效率是关键。在计算单个距离后添加所有段每条路径需要 7 毫秒。我想让它至少快 1 毫秒。
编辑： 我需要使用 WGS84 椭球来计算距离，因此球面 (Haversine) 是不够的。我想我可以以 1m 的精度工作。点沿路径随机分布。有的可以20公里，有的不到1公里。
这是带点的路径（十进制lat/lon）：
[(49.009722, 2.547778), (49.015556, 2.573611), (49.021389, 2.599167), (49.039167, 2.676389), (49.048056, 2.715), (49.044444, 2.835), (49.041667, 2.928333), (49.042778, 2.942222), (49.051667, 3.066667), (49.061389, 3.205), (49.072222, 3.357222), (49.085, 3.536944), (49.086111, 3.550833), (49.097778, 3.729444), (49.113056, 3.963056), (49.130833, 4.238056), (49.138889, 4.361667), (49.1925, 4.564444), (49.306667, 4.995556), (49.333611, 5.096944), (49.395, 5.329167), (49.490556, 5.690833), (49.514444, 5.781111), (49.53, 5.845833), (49.599444, 6.127778), (49.637222, 6.281667), (49.673333, 6.440278), (50.0475, 8.078333), (50.053611, 8.637222), (50.056667, 8.800278), (50.063056, 9.19), (50.066944, 9.486389), (50.07, 9.783056), (50.072778, 10.098611), (50.073333, 10.242778), (50.075278, 10.728889), (50.046667, 10.863333), (50.0325, 10.930278), (49.981111, 11.172222), (49.969722, 11.225833), (49.961111, 11.491389), (49.959444, 11.547222), (49.957222, 11.617222), (49.946111, 11.9325), (49.937222, 12.343889), (49.933333, 12.47), (49.9325, 12.498056), (49.928611, 12.624167), (49.924167, 12.764444), (49.919444, 12.918611), (49.910833, 13.199167), (49.909444, 13.241111), (49.907778, 13.283611), (49.900556, 13.481944), (50.077222, 13.840278), (50.124167, 13.995556), (50.182778, 14.189722), (50.220278, 14.315), (50.268889, 14.478056), (50.211389, 14.403611), (50.166389, 14.345), (50.133611, 14.3025), (50.100833, 14.26)]
我发现了cartopy - a package providing geodesic calculations that uses shapely and proj。然而，关于 cartpopy 的文档没有给出相关示例，我被困在了这一点上。基本上 geometry_length 一次性给出了一个有形物体的长度，所以我用下面的方式来做：
#defining the geoid on which to make calculations
myGeod = geodesic.Geodesic(6378137.0,1 / 298.257223563)
#making my list of latlon (in decimal degrees) into a shapely
shapelyObject = LineString(list(latlon_dd))
#applying the method on the shapelyObject given the defined ellipsoid
myGeod.geometry_length(shapelyObject)
我想计算长度（以米为单位），应该在 917,315.3 米左右。相反，我得到这个 ValueError:
ValueError Traceback (most recent call last)
<ipython-input-243-7c75042775e3> in <module>
6
7 #applying the method on the shapelyObject given the defined ellipsoid
----> 8 myGeod.geometry_length(shapelyObject)
lib\cartopy\geodesic\_geodesic.pyx in cartopy.geodesic._geodesic.Geodesic.geometry_length()
lib\cartopy\geodesic\_geodesic.pyx in cartopy.geodesic._geodesic.Geodesic.geometry_length()
lib\cartopy\geodesic\_geodesic.pyx in cartopy.geodesic._geodesic.Geodesic.inverse()
ValueError: Expecting input points to be (N, 2), got (1, 63)
提前致谢！
将根据找到的回复回答我的问题 here。显然这是一个错误，目前的解决方法是使用：
myGeod.geometry_length(np.array(shapelyObject.coords))
而不是
myGeod.geometry_length(shapelyObject)
将在最终解决方案可用时更新。
您可能不想包含另一个库，但是 LatLon 有一个内置方法来计算两个 lat/lon 点之间的 WGS84 距离（以千米为单位）。
示例取自链接页面：
>> palmyra = LatLon(Latitude(5.8833), Longitude(-162.0833)) # Location of Palmyra Atoll
>> honolulu = LatLon(Latitude(21.3), Longitude(-157.8167)) # Location of Honolulu, HI
>> distance = palmyra.distance(honolulu) # WGS84 distance in km
>> print distance
1766.69130376
>> print palmyra.distance(honolulu, ellipse = 'sphere') # FAI distance in km
1774.77188181
编辑：发帖后才注意到你想计算一条路径的距离，其中有很多点不在两个点之间（一次），抱歉...
您仍然可以使用 pyproj 并一次解决整个列表
geod = pyproj.Geod(ellps='WGS84')
_, _, distances_in_meters = geod.inv(
lons1_float_or_list_or_numpy_array,
lats1_float_or_list_or_numpy_array,
lons2_float_or_list_or_numpy_array,
lats2_float_or_list_or_numpy_array)
要从您提供的内容中获取此格式，我们可以使用 4 个理解和 ittertools.combinations。
your_list = [(49.009722, 2.547778), (49.015556, 2.573611), ...]
edges = ittertools.combinations(your_list, 2)
lons1 = [edge[0][1], for edge in edges]
lats1 = [edge[0][0], for edge in edges]
lons2 = [edge[1][1], for edge in edges]
lats2 = [edge[1][0], for edge in edges]