Haversine公式（半正矢公式）

原理分析

直接套公式，簡潔高效，

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

Python代碼

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 ∣ ) \cos\theta = \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

Python代碼

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度，所以在赤道上經度相差一度為40075020/360=111319.5m=111.3195km；對于緯度不為0的，經度相差一度為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 & =r\sin\theta \cos\phi \\ y & =r\sin\theta \sin\phi \\ z & =r\cos\theta \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}{r\sin\theta})=\arcsin(\frac{y}{r\sin\theta})=\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，2\pi] 徑向距離r∈[0，+∞]， 傾角(天頂角)θ∈[0，π]， 方位角?∈[0，2π]