快速半正弦逼近（Python/Pandas）

Posted 2023-02-16

【中文标题】快速半正弦逼近（Python/Pandas）

【英文标题】：Fast Haversine Approximation (Python/Pandas)

【发布时间】：2015-06-15 05:51:38

【问题描述】：

Pandas 数据框中的每一行都包含 2 个点的 lat/lng 坐标。使用下面的 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'])

【问题讨论】：

比近似更好的方法是分析函数以准确了解为什么它需要太长时间，然后使用 ctypes/Cython/numba 将函数原样转换为 C 函数运行没有太多的开销。您可能需要修改调用约定以使用每个 pandas Series 数据列底层数据的 numpy 数组 values，您还可以签出 numpy.ctypeslib 以轻松从 numpy 数组转换为与 ctypes 兼容的数组。看起来很多，但实际上这是在 Python 中访问 C 函数的一种非常简单的方法。

可以避免对大多数候选人进行计算。计算距起点 50 英里的最小和最大经度和纬度。然后使用这些最小值和最大值来淘汰大多数候选人。

您还可以考虑从数据中构建一个 k-d 树，而不是将其存储在像 DataFrame 这样的关系结构中。那么获得给定点的邻居会很便宜，也许你只能按需计算距离。应用程序是否总是需要每一对？另一种选择可能是将点聚类并使用每个聚类的质心/平均值作为代理。然后任何两点之间的距离将近似为仅聚类中心之间的距离。不过，像这样的奇思妙想是否真的比蛮力更好，这是推测性的。

@Nyxynyx 您在问题中提供的函数给出了很大的圆距。您评论中的计算给出了欧几里得距离。因为地球的半径很大，小距离绝对可以用欧几里德的版本来近似。

是的，欧几里得近似适用于足够小的距离。你甚至不需要为此做一个apply，可以直接使用数据框中的列。

【参考方案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 = 6367 * 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 提供了对整个数据数组进行操作的函数，可以避免循环并显着提高性能。

这是vectorization 的示例。

【讨论】：

很高兴知道array programming这个术语，在MATLAB中没有遇到过。

非常感谢您。小建议：添加一个带有实际坐标而不是随机值的真实示例用法，以阐明输入格式。

请注意，当一对参数是 Series 而另一个是元组时，这也有效：haversine_np(pd.Series([-74.00594, -122.41942]), pd.Series([40.71278, 37.77493]), -87.65005, 41.85003) 计算（纽约、旧金山）和芝加哥之间的距离。

另一个小建议：您可能希望将函数参数的顺序交换为lat, lon。在许多来源中，纬度是第一位的，例如在en.wikipedia.org/wiki/Horizontal_position_representation.

我向 sklearn 提出了功能请求以添加您的代码：github.com/scikit-learn/scikit-learn/issues/17212

【参考方案2】：

纯粹为了一个说明性的例子，我在@ballsdotballs 的答案中采用了numpy 版本，并且还制作了一个配套的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\n", res);

    int i;
    for (i=0; i < arr_size; i++)
        printf("%3.3f, ", kms[i]);
    
    printf("\n");

请注意，我们正在努力遵守 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 之间随机选择 lats 和 lons（以度为单位），但这个解释并不重要。

接下来我们需要做的是编译我们的 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 加载的函数，就像它们是指针一样。代码如下：

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

将它们放在haversine.py 的__main__ 块中，整个文件现在看起来像这样：

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

要运行它，它将分别运行和计时 Python 和 ctypes 版本并打印一些结果，我们可以这样做

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 更好的选择。这正是 numpy 旨在解决的问题，因此，只要 (a) 在应用程序中合并 numpy 数据类型和 (b) 存在一种简单的方法，就可以自制您自己的 ctypes 代码将您的代码映射到 numpy 等效项，效率不高。

但是，当您更喜欢用 C 编写东西但在 Python 中调用它时，或者在不实际依赖 numpy 的情况下（在嵌入式系统中 @例如不能安装987654366@）。

【讨论】：

这太棒了！

【参考方案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))

【讨论】：

请注意参数顺序应该是纬度，经度与许多 GIS 库不同

【参考方案4】：

@derricw's vectorised solution 的一个微不足道的扩展，您可以使用 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】：

矢量化函数指定“所有参数必须等长”。根据this，通过扩展“更大”数据集的边界，可以有效地找到所有 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】：

其中一些答案是“四舍五入”地球的半径。如果您将这些与其他距离计算器（例如 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))