根据经纬度计算两地之间的距离

• Haversine公式（半正矢公式）
• • 原理分析
  • Python代码
• 向量法求两地距离
• • 原理分析
  • Python代码
• 知识补充
• 球坐标系与直角坐标系的转换

直接套公式，简洁高效。

推导过程：
https://zhuanlan.zhihu.com/p/202993882?utm_source=wechat_session
https://www.douban.com/note/766524321/

import math

R = 6371.393
Pi = math.pi


# A地
weiduA, jingduA = 39.890115, 116.295794

# B地
weiduB, jingduB = 31.212544, 121.400053


a = (math.sin(math.radians(weiduA/2-weiduB/2)))**2
b = math.cos(weiduA*Pi/180) * math.cos(weiduB*Pi/180) * (math.sin((jingduA/2-jingduB/2)*Pi/180))**2

L = 2 * R * math.asin((a+b)**0.5)

print(L)

思路：以地心为原点，建立空间直角坐标系。将地球上任意两点的经纬度坐标转为直角坐标。然后用向量求出夹角余弦值，最后用角度求弧长。
误差：地球不是完美的球体。

以地心为原点，建立空间直角坐标系：

X轴为地心到（经度：0°、纬度：0°）的向量
Y轴为地心到（经度：90°、纬度：0°）的向量
Z轴为地心到（纬度：90°）的向量

由此可以得到由经纬度到空间直角坐标系的对应关系：

x = cos ⁡ ( 纬 度 ) cos ⁡ ( 经 度 )   ， y = cos ⁡ ( 纬 度 ) sin ⁡ ( 经 度 )   ， z = sin ⁡ ( 纬 度 ) x = cos(纬度)cos(经度) ，y = cos(纬度) sin(经度) ，z = sin(纬度) x=cos(纬度)cos(经度) ，y=cos(纬度)sin(经度) ，z=sin(纬度)

由两个点对应的向量求出向量夹角：
cos ⁡ θ = V A ⋅ V B ∣ V A ∣ × ∣ V B ∣ θ = arccos ⁡ ( V A ⋅ V B ∣ V A ∣ ⋅ ∣ V B ∣ ) costheta = frac{V_A cdot V_B }{|V_A| times |V_B|} qquad qquad theta=arccos(frac{V_A cdot V_B }{|V_A| cdot |V_B|}) cosθ=∣VA​∣×∣VB​∣VA​⋅VB​​θ=arccos(∣VA​∣⋅∣VB​∣VA​⋅VB​​)

假设地球为理想球体：半径大约3959英里(6371.393千米) 。这个数字是地心到地球表面所有各点距离的平均值，平均半径=(赤道半径×2+极半径)/3。

则由半径和夹角可求弧长：
弧 长 l = π R θ 180 = α R 弧长l=pi R frac{theta}{180} = alpha R 弧长l=πR180θ​=αR

其中 θ theta θ是圆心角度数（角度制），R是半径，L是圆心角弧长， α alpha α是圆心角度数（弧度制）。 α = π θ / 180 alpha = pi theta/180 α=πθ/180

import math

R = 6371.393
Pi = math.pi


# A地
weiduA, jingduA = 39.890115, 116.295794
# 转为空间直角坐标
xA = math.cos(math.radians(weiduA))*math.cos(math.radians(jingduA))
yA = math.cos(math.radians(weiduA))*math.sin(math.radians(jingduA))
zA = math.sin(math.radians(weiduA))


# B地
weiduB, jingduB = 31.212544, 121.400053
# 转为空间直角坐标
xB = math.cos(weiduB*Pi/180) * math.cos(jingduB*Pi/180)
yB = math.cos(weiduB*Pi/180) * math.sin(jingduB*Pi/180)
zB = math.sin(weiduB*Pi/180)


# 开始计算
cosalpha = (xA*xB+yA*yB+zA*zB)/((xA*xA+yA*yA+zA*zA)*(xB*xB+yB*yB+zB*zB))**0.5

alpha = math.acos(cosalpha)

L = alpha * R

print(L)

• 在经线上纬度差1度对应的实际距离是111.2018千米
• 在赤道上经度差1度对应的实际距离是111.3195千米
• 在除赤道外的其他纬线上，经度差1度对应的实际距离是111.3195*cos纬度

从理论上讲，全部的经线长度都相等，无论沿哪条经线南北极之间的距离都相等。
所以从理论上算，一条经线的长度=平均半径乘以圆周率= π pi πR=20016.321441933433千米。所以，在同一条经线上，纬度差1度对应的实际距离是 π pi πR/180=111.20178578851908千米。

赤道周长: 40075020m。
因为赤道被分为了360度，所以在赤道上经度差1度对应的实际距离是40075020/360=111319.5m=111.3195km；对于纬度不为0的情况，在同一条纬线上，经度差1度对应的实际距离是111.3195 × times ×cos纬度。
（另一种理论计算结果应该是111.2018 × times ×cos纬度，因为平均半径R × times ×cos纬度等于该纬度对应的小圆半径（一条纬线就是一个圈，小圆说的就是这个圈），1度所对应的弧长就是2 π pi πR × times ×cos纬度/360=111.2018 × times ×cos纬度）

总结，对于日常的学习生活来说，相差一度取111km、圆周率用3.14就够了，不必太过较真。

x = r sin ⁡ θ cos ⁡ ϕ y = r sin ⁡ θ sin ⁡ ϕ z = r cos ⁡ θ ⟺ r = x 2 + y 2 + z 2 θ = arccos ⁡ ( z r ) = arcsin ⁡ ( x 2 + y 2 r ) = arctan ⁡ ( x 2 + y 2 z ) ϕ = arccos ⁡ ( x r sin ⁡ θ ) = arcsin ⁡ ( y r sin ⁡ θ ) = arctan ⁡ ( y x ) begin{aligned} begin{aligned} x & =rsintheta cosphi \ y & =rsintheta sinphi \ z & =rcostheta end{aligned} qquad Longleftrightarrow qquad begin{aligned} r & =sqrt{x^2+y^2+z^2} \ theta & =arccos(frac{z}{r})=arcsin(frac{sqrt{x^2+y^2}}{r})=arctan(frac{sqrt{x^2+y^2}}{z}) \ phi & =arccos(frac{x}{rsintheta})=arcsin(frac{y}{rsintheta})=arctan(frac{y}{x}) end{aligned} end{aligned} xyz​=rsinθcosϕ=rsinθsinϕ=rcosθ​⟺rθϕ​=x2+y2+z2 ​=arccos(rz​)=arcsin(rx2+y2 ​​)=arctan(zx2+y2 ​​)=arccos(rsinθx​)=arcsin(rsinθy​)=arctan(xy​)​​

径 向 距 离 r ∈ [ 0 ， + ∞ ] ，   倾 角 ( 天 顶 角 ) θ ∈ [ 0 ， π ] ，   方 位 角 ϕ ∈ [ 0 ， 2 π ] 径向距离r in [0，+infty ]， 倾角(天顶角)theta in [0，pi ]， 方位角phi in [0，2pi] 径向距离r∈[0，+∞]， 倾角(天顶角)θ∈[0，π]， 方位角ϕ∈[0，2π]