在Python中利用距离矩阵计算经纬点间的距离

import math

def haversine(lat1, lon1, lat2, lon2):
    R = 6371  # 地球平均半径，单位为公里
    dlat = math.radians(lat2 - lat1)
    dlon = math.radians(lon2 - lon1)
    a = math.sin(dlat / 2) * math.sin(dlat / 2) + math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) * math.sin(dlon / 2) * math.sin(dlon / 2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    distance = R * c
    return distance

import math

def vincenty(lat1, lon1, lat2, lon2):
    a = 6378137  # 赤道半径，单位为米
    b = 6356752.3142  # 极半径，单位为米
    f = 1 / 298.257223563  # 扁率
    L = math.radians(lon2 - lon1)
    U1 = math.atan((1 - f) * math.tan(math.radians(lat1)))
    U2 = math.atan((1 - f) * math.tan(math.radians(lat2)))
    sinU1 = math.sin(U1)
    cosU1 = math.cos(U1)
    sinU2 = math.sin(U2)
    cosU2 = math.cos(U2)
    lambdaL = L
    lambdaP = 2 * math.pi
    iterLimit = 100
    while abs(lambdaL - lambdaP) > 1e-12 and iterLimit > 0:
        sinLambdaL = math.sin(lambdaL)
        cosLambdaL = math.cos(lambdaL)
        sinSigma = math.sqrt((cosU2 * sinLambdaL) ** 2 + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambdaL) ** 2)
        cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambdaL
        sigma = math.atan2(sinSigma, cosSigma)
        sinAlpha = cosU1 * cosU2 * sinLambdaL / sinSigma
        cosSqAlpha = 1 - sinAlpha ** 2
        cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha
        C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
        lambdaP = lambdaL
        lambdaL = L + (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM ** 2)))
        iterLimit -= 1
    if iterLimit == 0:
        return float('inf')  # 超过迭代次数，认为计算失败
    uSq = cosSqAlpha * (a ** 2 - b ** 2) / (b ** 2)
    A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
    B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
    deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM ** 2) - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma ** 2) * (-3 + 4 * cos2SigmaM ** 2)))
    distance = b * A * (sigma - deltaSigma)
    return distance

distance = haversine(39.9075, 116.39723, 31.2304, 121.4737)
print(distance)