怎么用Python计算地球上任意两个用经纬度表示的点的弧面距离？

这是来自知乎上的问题。问我的时候，恰好我在写一个和向量计算相关的文章，于是灵光乍现，顺手写了这样一个答案。该算法未经严格验证，请谨慎参考。具体思路如下。

将两个点的经纬度换算成空间坐标；计算地心与两个点所成的两个向量的点积；点积除以两个向量的模（也就是地球半径）之积，结果就是向量夹角的余弦；反余弦值对应着两点所在大圆（即经过两点的地球表面最大的圆）的弧度；弧度乘以地球半径，即得弧长。

代码如下。

>>> def get_arc(p0, p1, r=1):z0 = r*np.sin(np.radians(p0[1]))x0 = r*np.cos(np.radians(p0[1]))*np.cos(np.radians(p0[0]))y0 = r*np.cos(np.radians(p0[1]))*np.sin(np.radians(p0[0]))z1 = r*np.sin(np.radians(p1[1]))x1 = r*np.cos(np.radians(p1[1]))*np.cos(np.radians(p1[0]))y1 = r*np.cos(np.radians(p1[1]))*np.sin(np.radians(p1[0]))theta = np.arccos(np.dot((x0,y0,z0),(x1,y1,z1))/(r*r))return theta * r>>> r = 6377.830 # 使用赤道半径，单位：km>>> 北京 = (116.5,40.0)>>> 济南 = (117.0,36.4)>>> 悉尼 = (151.2,-33.9)>>> 纽约 = (-74.0,40.5)>>> 巴西利亚 = (-48.0,-15.9)>>> get_arc(北京, 济南, r=r)403.10853123743505>>> get_arc(北京, 悉尼, r=r)8966.0798138>>> get_arc(北京, 纽约, r=r)11012.732277406261>>> get_arc(北京, 巴西利亚, r=r)16962.014473669606