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