下面是同一函数的向量化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

输入都是值的数组，它应该能够立即处理数百万个点。要求是输入是ndarrays，但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的例子。

为了便于说明，我在@ballsdotballs的答案中采用了numpy版本，并通过ctypes调用了一个配套的C实现。由于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和隆恩（以度为单位），但这对这个解释来说并不重要。
我们需要做的下一件事是编译我们的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）存在一种简单的方法将代码Map到numpy等价物，自制自己的ctypes代码都不是很有效。
但是，当您喜欢用C编写但用Python调用某些内容时，或者依赖numpy不切实际的情况下（例如，在无法安装numpy的嵌入式系统中），了解如何做到这一点仍然非常有帮助。

如果允许使用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))

对@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

vectorized函数指定“所有参数必须长度相等”。通过扩展“较大”数据集的边界，根据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))

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

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