经纬度计算距离公式(测地线距离)

前言

做的几个项目都和两个poi测地线距离相关，特地总结一下，方便使用。

计算测地线距离

求上海和西安的距离。

# (纬度,经度)
point1 = (31.231706, 121.472644) # 上海
point2 = (34.263161, 108.948024) # 西安

方法1: 调用geopy库函数

geopy 可以使用第三方地理编码器和其他数据源定位全球地址、城市、国家和地标的坐标。它的一个重要功能就是实现经纬度地理位置转换，即求两点之间的测地线距离。但geopy速度相对较慢。

# geopy版本为'2.2.0'
import geopy.distance
# 输出距离,单位为km
geopy.distance.geodesic(point1, point2).km
"""
输出:
1220.029989895109
"""

方法2: 使用半正矢公式

半正矢公式是一种根据两点的经度和纬度来确定球面两点之间距离的计算方法。原理参见 wiki百科-Haversine formula。

from math import radians, sin, cos, asin, sqrt
def get_geodesic(point1, point2):
    """
    point1,point2格式为(纬度,经度)
    如(121.472644, 31.231706)
    """
    lat1,lon1 = map(radians, point1)
    lat2,lon2 = map(radians, point2)
    # 半正矢公式(经纬度求解两点球面距离)
    # Haversine formula
    diff_lon = lon2 - lon1 
    diff_lat = lat2 - lat1 
    P = sin(diff_lat/2)**2 + cos(lat1)*cos(lat2)*sin(diff_lon/2)**2
    Q = 2 * asin(sqrt(P))   
    # 地球半径(km)
    R = 6371  
    # 测地线距离
    geodesic = Q * R
    return geodesic
get_geodesic(point1, point2)
"""
输出:
1217.912379941211
"""

注意，方法1和方法2之间存在一定误差，大量的数据表明，半正矢公式相比于geopy.distance.geodesic的平均误差为 0.16%，即1000km相差约1.6km，见下一节。

geopy和半正矢公式误差比较

平均误差计算：

import numpy as np
from tqdm import tqdm
L = []
for i in tqdm(range(10000)):
    # 纬度范围为[-90,90],经度范围为[-180,180]
    lat1, lat2 = np.random.random(2)*np.sign(np.random.randn(2))*90
    lon1, lon2 = np.random.random(2)*np.sign(np.random.randn(2))*180
    point1 = (lat1, lon1)
    point2 = (lat2, lon2)
    d1 = geopy.distance.geodesic(point1, point2).km
    d2 = get_geodesic(point1, point2)
    L.append((d1, d2))
L1 = [abs(x[0]-x[1])/x[0] for x in L]
print(np.mean(L1))
"""
输出:
0.0016
"""

geopy和半正矢公式耗时比较

在AMD Ryzen 7 5800U上测试了两个函数的耗时，其结果如下：

geopy.distance.geodesic(point1, point2).km 万次计算平均每次耗时

155 µs ± 3.97 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

get_geodesic(point1, point2) 百万次计算平均每次耗时

1.08 µs ± 14.4 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

可见，geopy时间复杂度较大，单次运行时长是半正矢公式的143.5倍。

附

计算两个坐标之间的距离资料：

• https://www.hhlink.com/%E7%BB%8F%E7%BA%AC%E5%BA%A6/
• 谷歌地球