如何使用带有 NumPy 数组的矢量化来使用 Geopy 库为大型数据集计算测地线距离？答案

【问题标题】：How to use Vectorization with NumPy arrays to calculate geodesic distance using Geopy library for a large dataset?如何使用带有 NumPy 数组的矢量化来使用 Geopy 库为大型数据集计算测地线距离？
【发布时间】：2018-10-20 20:26:39
【问题描述】：

我正在尝试从一个数据帧计算测地距离，该数据帧由四列纬度和经度数据组成，大约有 300 万行。我使用 apply lambda 方法来完成，但完成任务需要 18 分钟。有没有办法使用带有 NumPy 数组的矢量化来加快计算速度？谢谢你的回答。

我的代码使用 apply 和 lambda 方法：

from geopy import distance

df['geo_dist'] = df.apply(lambda x: distance.distance(
                              (x['start_latitude'], x['start_longitude']),
                              (x['end_latitude'], x['end_longitude'])).miles, axis=1)

更新：

我正在尝试这段代码，但它给了我错误：ValueError：具有多个元素的数组的真值不明确。使用 a.any() 或 a.all()。如果有人可以提供帮助，请不胜感激。

df['geo_dist'] = distance.distance(
                          (df['start_latitude'].values, df['start_longitude'].values),
                          (df['end_latitude'].values, df['end_longitude'].values)).miles

【问题讨论】：

Geopy 距离例程目前不支持矢量化。也许测地线例程的作者@cffk 可能会在那里提出解决方案？此问题正在github.com/geopy/geopy/issues/189 中进行跟踪

标签： python numpy coordinates vectorization geopy

【解决方案1】：

我认为您可能会考虑为此使用 geopandas，它是 pandas（因此是 numpy）的扩展，旨在非常快速地执行这些类型的计算。

具体来说是has a method for calculating the distance between sets of points in a GeoSeries，可以是GeoDataFrame的一列。我相当肯定这种方法利用numexpr 进行矢量化。

它应该看起来像这样，您将数据框转换为具有（至少）两个 GeoSeries 列的 GeoDataFrame，您可以将它们用于起点和终点。这应该返回一个GeoSeries 对象：

import pandas as pd
import geopandas as gpd
from shapely.geometry import Point

geometry = [Point(xy) for xy in zip(df.longitude, df.latitude)]
gdf = gpd.GeoDataFrame(df, crs={'init': 'epsg:4326'}, geometry=geometry)

distances = gdf.geometry.distance(gdf.destination_geometry)

【讨论】：

gdf.destination_geometry 从何而来？我通过使用一组不同的经度和纬度创建 dest_geometry，创建 dest_gdf，然后调用 gdf.geometry.distance(dest_gdf.geometry) 来使其工作。此外，该 gpd.GeoDataFrame 行在授权码上给出了 FutureWarning。尽我所能告诉该行应该是： gdf = gpd.GeoDataFrame(df, crs='epsg:4326', geometry=geometry)

【解决方案2】：

用 numpy 来回走动：

from geopy import distance

lats = df['latitude'].values
lons = df['longitude'].values
latsNext = np.roll(lats, 1)
lonsNext = np.roll(lons, 1)
dists = [distance.distance((lat0, lon0),(lat1, lon1)).kilometers for lat0, lon0, lat1, lon1 in zip(lats, lons, latsNext, lonsNext)]
dists = np.roll(dists, -1)
dists[-1] = np.nan
df['distance'] = dists

【讨论】：

【解决方案3】：

你的问题的答案：你不能用geopy做你想做的事。我不熟悉这个包，但是错误回溯表明这个函数以及这个包中可能的所有其他函数都没有考虑到矢量化计算。

现在，如果您可以处理大圆距离，那么我建议您尝试使用 astropy.coordinates 包，我能够以矢量方式计算点之间的 separations。

这是一个基于我对另一个问题的回答的示例：Finding closest point:

from astropy.units import Quantity
from astropy.coordinates import SkyCoord, EarthLocation
from astropy.constants import R_earth
import numpy as np

lon1 = Quantity([-71.312796, -87.645307, -87.640426, -87.635513,
                 -87.630629, -87.625793 ], unit='deg')
lat1 = Quantity([41.49008, 41.894577, 41.894647, 41.894713,
                 41.894768, 41.894830], unit='deg')
lon2 = Quantity([-81.695391, -87.645307 + 0.5, -87.640426, -87.635513 - 0.5,
                 -87.630629 + 1.0, -87.625793 - 1.0], unit='deg')
lat2 = Quantity([41.499498, 41.894577 - 0.5, 41.894647, 41.894713 - 0.5,
                 41.894768 - 1.0, 41.894830 + 1.0], unit='deg')

pts1 = SkyCoord(EarthLocation.from_geodetic(lon1, lat1, height=R_earth).itrs, frame='itrs')
pts2 = SkyCoord(EarthLocation.from_geodetic(lon2, lat2, height=R_earth).itrs, frame='itrs')

然后，两组点之间的距离可以计算为：

>>> dist = pts2.separation(pts1)
>>> print(dist)
<Angle [ 7.78350849, 0.62435354, 0., 0.62435308, 1.25039805, 1.24353876] deg>

近似转换为距离：

>>> np.deg2rad(pts2.separation(pts1)) * R_earth / u.rad
<Quantity [ 866451.17527216,  69502.31527953,      0.        ,
             69502.26348614, 139192.86680148, 138429.29874024] m>

将第一个值与您从geopy 的示例中得到的值进行比较：

>>> distance.distance((41.49008, -71.312796), (41.499498, -81.695391)).meters
866455.4329098687

编辑：实际上，这很可能实际上会为您提供您所追求的测地距离，但请务必检查the description of EarthLocation。

【讨论】：