纯粹为了一个说明性的例子，我在@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@）。