根据经纬度计算地球距离-Numpy向量版

输入两个numpy数组,shape都是(n,2)，每一行代表一个点，例如：

latlon1 = array([[ 22.7373478, 113.8013962],
       [ 22.6623782, 113.8286608]])

latlon2 = array([[ 22.6623782, 113.8286608],
       [ 22.6874816, 113.8103144]])

注意：两个数组需要满足行数一样或者其中一个数组的行数为1，当行数一样时，计算的是两个数组同一行的点的距离；有一个数组行数是1时，计算的是该点与另一数组所有点的距离。
代码如下：

def haversine(latlon1, latlon2):
    """
    计算两个经纬度点向量之间的地球距离
    输入：向量1（纬度，经度），向量2（纬度，经度）
    输出：距离向量
    """
    latlon1 = latlon1.reshape(-1,2) #固定形状
    latlon2 = latlon2.reshape(-1,2) #固定形状
    lat1, lon1 = latlon1[:,1],latlon1[:,0]
    lat2, lon2 = latlon2[:,1],latlon2[:,0]
    
    lon1 = list(map(radians,lon1))
    lat1 = list(map(radians,lat1))
    lon2 = list(map(radians,lon2))
    lat2 = list(map(radians,lat2))

    dlon = np.array(lon2) - np.array(lon1)
    dlat = np.array(lat2) - np.array(lat1)

    sindlat_p2 = list(map(lambda x:x**2,list(map(sin,(dlat/2).tolist()))))
    coslat1 = list(map(cos,lat1))
    coslat2 = list(map(cos,lat2))
    sindlon_p2 = list(map(lambda x:x**2,list(map(sin,(dlon/2).tolist()))))
    a = np.array(sindlat_p2)+np.array(coslat1)*np.array(coslat2)*np.array(sindlon_p2)
    c = np.array(list(map(asin,list(map(sqrt,a.tolist())))))*2
    r = 6370996.81  # 地球半径
    distance = c*r
    
    return distance

以上数据的计算结果：

array([8792.90409277, 3366.74435249])