快速半正弦近似（Python/Pandas）
- 2025-01-09 08:46:00
- admin 原创
- 11
问题描述：
Pandas 数据框中的每一行包含 2 个点的纬度/经度坐标。使用下面的 Python 代码，计算许多（数百万）行的这两个点之间的距离需要很长时间！
考虑到 2 个点相距不到 50 英里且精度不是很重要，是否有可能加快计算速度？
from math import radians, cos, sin, asin, sqrt
def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a))
km = 6367 * c
return km
for index, row in df.iterrows():
df.loc[index, 'distance'] = haversine(row['a_longitude'], row['a_latitude'], row['b_longitude'], row['b_latitude'])
解决方案 1：
以下是同一函数的矢量化 numpy 版本：
import numpy as np
def haversine_np(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
All args must be of equal length.
"""
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2
c = 2 * np.arcsin(np.sqrt(a))
km = 6378.137 * c
return km
输入都是值数组，它应该能够立即处理数百万个点。要求是输入是 ndarray，但 pandas 表的列可以正常工作。
例如，使用随机生成的值：
>>> import numpy as np
>>> import pandas
>>> lon1, lon2, lat1, lat2 = np.random.randn(4, 1000000)
>>> df = pandas.DataFrame(data={'lon1':lon1,'lon2':lon2,'lat1':lat1,'lat2':lat2})
>>> km = haversine_np(df['lon1'],df['lat1'],df['lon2'],df['lat2'])
或者如果您想创建另一列：
>>> df['distance'] = haversine_np(df['lon1'],df['lat1'],df['lon2'],df['lat2'])
在 Python 中循环遍历数据数组非常慢。Numpy 提供了对整个数据数组进行操作的函数，让您可以避免循环并大幅提高性能。
这是矢量化的一个例子。
解决方案 2：
纯粹为了举例说明，我采用了numpy@ballsdotballs 答案中的版本，同时还制作了一个配套的 C 实现，通过 进行调用ctypes。由于numpy是一种高度优化的工具，我的 C 代码不太可能如此高效，但应该比较接近。这里最大的优势是，通过运行一个包含 C 类型的示例，它可以帮助您了解如何将您自己的个人 C 函数连接到 Python，而开销又不会太大。当您只想通过在 C 源代码而不是 Python 中编写一小段来优化较大计算的一小部分时，这尤其好。numpy大多数情况下，简单地使用 就能解决问题，但对于那些您并不真正需要全部numpy并且不想添加耦合以要求numpy在整个代码中使用数据类型的情况，知道如何使用内置ctypes库并自己动手会非常方便。
首先让我们创建 C 源文件，名为haversine.c：
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
int haversine(size_t n,
double *lon1,
double *lat1,
double *lon2,
double *lat2,
double *kms){
if ( lon1 == NULL
|| lon2 == NULL
|| lat1 == NULL
|| lat2 == NULL
|| kms == NULL){
return -1;
}
double km, dlon, dlat;
double iter_lon1, iter_lon2, iter_lat1, iter_lat2;
double km_conversion = 2.0 * 6367.0;
double degrees2radians = 3.14159/180.0;
int i;
for(i=0; i < n; i++){
iter_lon1 = lon1[i] * degrees2radians;
iter_lat1 = lat1[i] * degrees2radians;
iter_lon2 = lon2[i] * degrees2radians;
iter_lat2 = lat2[i] * degrees2radians;
dlon = iter_lon2 - iter_lon1;
dlat = iter_lat2 - iter_lat1;
km = pow(sin(dlat/2.0), 2.0)
+ cos(iter_lat1) * cos(iter_lat2) * pow(sin(dlon/2.0), 2.0);
kms[i] = km_conversion * asin(sqrt(km));
}
return 0;
}
// main function for testing
int main(void) {
double lat1[2] = {16.8, 27.4};
double lon1[2] = {8.44, 1.23};
double lat2[2] = {33.5, 20.07};
double lon2[2] = {14.88, 3.05};
double kms[2] = {0.0, 0.0};
size_t arr_size = 2;
int res;
res = haversine(arr_size, lon1, lat1, lon2, lat2, kms);
printf("%d
", res);
int i;
for (i=0; i < arr_size; i++){
printf("%3.3f, ", kms[i]);
}
printf("
");
}
请注意，我们试图遵循 C 惯例。通过引用显式传递数据参数，使用size_t大小变量，并期望我们的haversine函数通过改变传递的输入之一来工作，以便在退出时包含预期的数据。该函数实际上返回一个整数，这是一个成功/失败标志，可供该函数的其他 C 级使用者使用。
我们需要找到一种方法来处理 Python 内部所有这些特定于 C 的小问题。
接下来，让我们将我们的numpy函数版本以及一些导入和一些测试数据放入一个名为的文件中haversine.py：
import time
import ctypes
import numpy as np
from math import radians, cos, sin, asin, sqrt
def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = (np.sin(dlat/2)**2
+ np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2)
c = 2 * np.arcsin(np.sqrt(a))
km = 6367 * c
return km
if __name__ == "__main__":
lat1 = 50.0 * np.random.rand(1000000)
lon1 = 50.0 * np.random.rand(1000000)
lat2 = 50.0 * np.random.rand(1000000)
lon2 = 50.0 * np.random.rand(1000000)
t0 = time.time()
r1 = haversine(lon1, lat1, lon2, lat2)
t1 = time.time()
print t1-t0, r1
我选择将纬度和经度（以度为单位）随机选择在 0 到 50 之间，但这对于这个解释来说并不重要。
接下来我们需要做的是编译我们的 C 模块，以便它可以被 Python 动态加载。我使用的是 Linux 系统（您可以在 Google 上轻松找到其他系统的示例），所以我的目标是编译haversine.c成共享对象，如下所示：
gcc -shared -o haversine.so -fPIC haversine.c -lm
我们还可以编译为可执行文件并运行它来查看C程序的main函数显示的内容：
> gcc haversine.c -o haversine -lm
> ./haversine
0
1964.322, 835.278,
现在我们已经编译了共享对象haversine.so，我们可以用ctypes它在 Python 中加载它，并且我们需要提供文件的路径来执行此操作：
lib_path = "/path/to/haversine.so" # Obviously use your real path here.
haversine_lib = ctypes.CDLL(lib_path)
现在haversine_lib.haversine它的作用与 Python 函数非常相似，只是我们可能需要进行一些手动类型编组以确保正确解释输入和输出。
numpy实际上，它提供了一些很好的工具，我将在这里使用numpy.ctypeslib。我们将构建一个指针类型，允许我们将numpy.ndarrays这些ctypes-loaded 函数作为指针传递。以下是代码：
arr_1d_double = np.ctypeslib.ndpointer(dtype=np.double,
ndim=1,
flags='CONTIGUOUS')
haversine_lib.haversine.restype = ctypes.c_int
haversine_lib.haversine.argtypes = [ctypes.c_size_t,
arr_1d_double,
arr_1d_double,
arr_1d_double,
arr_1d_double,
arr_1d_double]
请注意，我们告诉haversine_lib.haversine函数代理根据我们想要的类型来解释其参数。
现在，为了从 Python 中测试它，剩下的就是创建一个大小变量和一个将被变异的数组（就像在 C 代码中一样）以包含结果数据，然后我们可以调用它：
size = len(lat1)
output = np.empty(size, dtype=np.double)
print "====="
print output
t2 = time.time()
res = haversine_lib.haversine(size, lon1, lat1, lon2, lat2, output)
t3 = time.time()
print t3 - t2, res
print type(output), output
把所有内容放在__main__块中haversine.py，整个文件现在如下所示：
import time
import ctypes
import numpy as np
from math import radians, cos, sin, asin, sqrt
def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = (np.sin(dlat/2)**2
+ np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2)
c = 2 * np.arcsin(np.sqrt(a))
km = 6367 * c
return km
if __name__ == "__main__":
lat1 = 50.0 * np.random.rand(1000000)
lon1 = 50.0 * np.random.rand(1000000)
lat2 = 50.0 * np.random.rand(1000000)
lon2 = 50.0 * np.random.rand(1000000)
t0 = time.time()
r1 = haversine(lon1, lat1, lon2, lat2)
t1 = time.time()
print t1-t0, r1
lib_path = "/home/ely/programming/python/numpy_ctypes/haversine.so"
haversine_lib = ctypes.CDLL(lib_path)
arr_1d_double = np.ctypeslib.ndpointer(dtype=np.double,
ndim=1,
flags='CONTIGUOUS')
haversine_lib.haversine.restype = ctypes.c_int
haversine_lib.haversine.argtypes = [ctypes.c_size_t,
arr_1d_double,
arr_1d_double,
arr_1d_double,
arr_1d_double,
arr_1d_double]
size = len(lat1)
output = np.empty(size, dtype=np.double)
print "====="
print output
t2 = time.time()
res = haversine_lib.haversine(size, lon1, lat1, lon2, lat2, output)
t3 = time.time()
print t3 - t2, res
print type(output), output
要运行它，它将ctypes分别运行和计时 Python 和版本并打印一些结果，我们只需这样做
python haversine.py
显示内容：
0.111340045929 [ 231.53695005 3042.84915093 169.5158946 ..., 1359.2656769
2686.87895954 3728.54788207]
=====
[ 6.92017600e-310 2.97780954e-316 2.97780954e-316 ...,
3.20676686e-001 1.31978329e-001 5.15819721e-001]
0.148446083069 0
<type 'numpy.ndarray'> [ 231.53675618 3042.84723579 169.51575588 ..., 1359.26453029
2686.87709456 3728.54493339]
正如预期的那样，该numpy版本稍微快一些（长度为 100 万的向量需要 0.11 秒），但我们的快速而粗糙的ctypes版本也毫不逊色：在相同数据上只需 0.148 秒。
让我们将其与 Python 中简单的 for 循环解决方案进行比较：
from math import radians, cos, sin, asin, sqrt
def slow_haversine(lon1, lat1, lon2, lat2):
n = len(lon1)
kms = np.empty(n, dtype=np.double)
for i in range(n):
lon1_v, lat1_v, lon2_v, lat2_v = map(
radians,
[lon1[i], lat1[i], lon2[i], lat2[i]]
)
dlon = lon2_v - lon1_v
dlat = lat2_v - lat1_v
a = (sin(dlat/2)**2
+ cos(lat1_v) * cos(lat2_v) * sin(dlon/2)**2)
c = 2 * asin(sqrt(a))
kms[i] = 6367 * c
return kms
当我将其放入与其他文件相同的 Python 文件中并对相同的百万元素数据进行计时时，我始终在我的计算机上看到大约 2.65 秒的时间。
因此，通过快速切换到，ctypes我们可以将速度提高约 18 倍。对于许多可以从访问裸的连续数据中受益的计算，您通常会看到比这更高的收益。
需要特别说明的是，我根本不认为这是比使用更好的选择。这正是要解决的numpy问题，因此，只要 (a) 在应用程序中包含数据类型是有意义的，并且 (b) 存在一种简单的方法将您的代码映射到等效代码中，那么编写自己的代码就不是很有效。numpy`ctypesnumpynumpy`
但是，当您喜欢用 C 编写某些东西但用 Python 调用它时，或者在依赖numpy不切实际的情况下（numpy例如，在无法安装的嵌入式系统中），知道如何做到这一点仍然非常有帮助。
解决方案 3：
如果允许使用 scikit-learn，我会给以下机会：
from sklearn.neighbors import DistanceMetric
dist = DistanceMetric.get_metric('haversine')
# example data
lat1, lon1 = 36.4256345, -5.1510261
lat2, lon2 = 40.4165, -3.7026
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
X = [[lat1, lon1],
[lat2, lon2]]
kms = 6367
print(kms * dist.pairwise(X))
解决方案 4：
这是@derricw 的矢量化解决方案的简单扩展，您可以使用它numba来将性能提高约 2 倍，而几乎无需更改代码。对于纯数值计算，这可能应该用于基准测试/测试，而不是可能更有效的解决方案。
from numba import njit
@njit
def haversine_nb(lon1, lat1, lon2, lat2):
lon1, lat1, lon2, lat2 = np.radians(lon1), np.radians(lat1), np.radians(lon2), np.radians(lat2)
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2
return 6367 * 2 * np.arcsin(np.sqrt(a))
与 Pandas 函数进行基准测试：
%timeit haversine_pd(df['lon1'], df['lat1'], df['lon2'], df['lat2'])
# 1 loop, best of 3: 1.81 s per loop
%timeit haversine_nb(df['lon1'].values, df['lat1'].values, df['lon2'].values, df['lat2'].values)
# 1 loop, best of 3: 921 ms per loop
完整基准测试代码：
import pandas as pd, numpy as np
from numba import njit
def haversine_pd(lon1, lat1, lon2, lat2):
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2
return 6367 * 2 * np.arcsin(np.sqrt(a))
@njit
def haversine_nb(lon1, lat1, lon2, lat2):
lon1, lat1, lon2, lat2 = np.radians(lon1), np.radians(lat1), np.radians(lon2), np.radians(lat2)
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2
return 6367 * 2 * np.arcsin(np.sqrt(a))
np.random.seed(0)
lon1, lon2, lat1, lat2 = np.random.randn(4, 10**7)
df = pd.DataFrame(data={'lon1':lon1,'lon2':lon2,'lat1':lat1,'lat2':lat2})
km = haversine_pd(df['lon1'], df['lat1'], df['lon2'], df['lat2'])
km_nb = haversine_nb(df['lon1'].values, df['lat1'].values, df['lon2'].values, df['lat2'].values)
assert np.isclose(km.values, km_nb).all()
%timeit haversine_pd(df['lon1'], df['lat1'], df['lon2'], df['lat2'])
# 1 loop, best of 3: 1.81 s per loop
%timeit haversine_nb(df['lon1'].values, df['lat1'].values, df['lon2'].values, df['lat2'].values)
# 1 loop, best of 3: 921 ms per loop
解决方案 5：
矢量化函数指定“所有参数必须长度相等”。通过扩展“较大”数据集的边界，根据此，可以有效地找到所有 i,j 对元素的距离。
from random import uniform
import numpy as np
def new_haversine_np(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
dlon = lon2 - lon1[:,None]
dlat = lat2 - lat1[:,None]
a = np.sin(dlat/2.0)**2 + np.cos(lat1[:,None]) * np.cos(lat2) * np.sin(dlon/2.0)**2
c = 2 * np.arcsin(np.sqrt(a))
km = 6367 * c
return km
lon1 = [uniform(-180,180) for n in range(6)]
lat1 = [uniform(-90, 90) for n in range(6)]
lon2 = [uniform(-180,180) for n in range(4)]
lat2 = [uniform(-90, 90) for n in range(4)]
new = new_haversine_np(lon1, lat1, lon2, lat2)
for i in range(6):
for j in range(4):
print(i,j,round(new[i,j],2))
解决方案 6：
Scikit-learn 库还有另一个用于计算半正矢距离的函数，称为haversine_distances函数，可用于查找两个坐标之间的距离，见以下示例：
from sklearn.metrics.pairwise import haversine_distances
import numpy as np
lat1, lon1 = [-34.83333, -58.5166646]
lat2, lon2 = [49.0083899664, 2.53844117956]
lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])
result = haversine_distances([[lat1, lon1], [lat2, lon2]])[0][1] # this formula returns a 2x2 array, this is the reason we used indexing.
print(result * 6373) # multiply by Earth radius to get kilometers
解决方案 7：
其中一些答案是“四舍五入”地球半径。如果你用其他距离计算器（如geopy）检查这些答案，这些函数将不起作用。
R=3959.87433如果您想要以英里为单位的答案， 则可以换成下面的转换常数。
如果您想要公里，请使用R= 6372.8。
lon1 = -103.548851
lat1 = 32.0004311
lon2 = -103.6041946
lat2 = 33.374939
def haversine(lat1, lon1, lat2, lon2):
R = 3959.87433 # this is in miles. For Earth radius in kilometers use 6372.8 km
dLat = radians(lat2 - lat1)
dLon = radians(lon2 - lon1)
lat1 = radians(lat1)
lat2 = radians(lat2)
a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2
c = 2*asin(sqrt(a))
return R * c
print(haversine(lat1, lon1, lat2, lon2))
相关推荐
热门文章
项目管理软件有哪些？
- 2024年20款好用的项目管理软件推荐，项目管理提效的20个工具和技巧
- 2024年开源项目管理软件有哪些？推荐5款好用的项目管理工具
- 项目管理软件有哪些？推荐7款超好用的项目管理工具
- 2024年常用的项目管理软件有哪些？推荐这10款国内外好用的项目管理工具
- 项目管理软件有哪些最好用？推荐6款好用的项目管理工具
- 项目管理软件哪个最好用？盘点推荐5款好用的项目管理工具
- 项目管理软件有哪些，盘点推荐国内外超好用的7款项目管理工具
- 项目管理软件排行榜：2024年项目经理必备5款开源项目管理软件汇总
- 2024项目管理软件排行榜（10类常用的项目管理工具全推荐）
- 项目管理必备:盘点2024年13款好用的项目管理软件
热门标签
云禅道AD
