python 计算经纬度距离，大地距离

python 计算经纬度距离，大地距离

最近做GPS轨迹挖掘，需要求两点间球面距离（大地距离）。网上没有很好的代码，于是写了这个。

先贴公式
对于点A 和 点B，
两点距离 D = arccos((sin北纬A×sin北纬B)＋(cos北纬A×cos北纬B×cosAB两地经度差))×地球平均半径 (Shormin)。 其中地球平均半径为6371.004 km，单位为km。

代码:

// An highlighted block
import numpy as np

def get_dis(pt1_:np.ndarray,pt2_:np.ndarray):
    """并行计算多点 距离"""
    R = 6371.004 # 地球半径
    pt1, pt2 = np.copy(pt1_[:,:2]), np.copy(pt2_[:,:2])
    pt1[:,1], pt2[:,1] = 90-pt1_[:,1], 90-pt2_[:,1] # 公式要求对纬度处理
    pt1, pt2 = np.radians(pt1), np.radians(pt2) # 角度转弧度
    a = np.sin(pt1[:,1]) * np.sin(pt2[:,1])* np.cos(pt1[:,0]-pt2[:,0])
    b = np.cos(pt1[:,1]) * np.cos(pt2[:,1])
    return R*np.arccos(a+b)

注意，是一次计算多点的哦。
比如

pt1 = np.array([[121,35], [121,36]])
pt2 = np.array([[122,36], [122,37]])
dis = get_dis(pt1, pt2) # km
print(dis) # [143.38274215 142.66524747]

顺便贴一个 连续轨迹 计算每点速度的代码

pts = np.array([[121,35,0], [121,36,2], [121,37,3]]) # [ [lon, lat, 时间戳]*n ]
pt1, pt2 = pts[:-1], pts[1:]
dis = get_dis(pt1,pt2)
time_delta = pt2[:,-1] - pt1[:,-1]
speed = dis/time_delta # km/s