根據經緯度計算兩地距離

最近工做須要，網上搜索了下根據經緯度計算兩地距離的方法，發現要麼是幾何法，畫圖、做一堆輔助線，而後證實推理，要麼二話不說直接套公式。這篇文章介紹一種容易理解的方式來求這個距離。

0b00 思路

地球是個不規則的橢球體、爲了簡便咱們看成球體來計算。
球體上兩地的最短距離就是通過兩點的大圓的劣弧長度。

思路以下:

弧長 ← 弦長(兩點距離) ← 兩點座標(直角座標) ← 經緯度

0b01 計算

1. 座標轉換

設

• 地球半徑爲 $R$
• 地心到 E 0° N 0° 的連線爲 x 軸
• 地心到 E 90° N 0° 的連線爲 y 軸
• 地心到 E 0° N 90° 的連線爲 z 軸
• 地球表面有一點 $A$, 經度爲 $e$, 緯度爲 $n$, 單位爲弧度

則 $A$ 的座標可表示爲:

$$ x = R \cdot cos(n) \cdot cos(e)\\ y = R \cdot cos(n) \cdot sin(e)\\ z = R \cdot sin(n) $$

代碼

const R = 6371
const {cos, sin, PI} = Math

let getPoint = (e, n) => {
    //首先將角度轉爲弧度
    e *= PI/180
    n *= PI/180
    reutrn {
        x: R*cos(n)*cos(e),
        y: R*cos(n)*sin(e),
        z: R*sin(n)
    }
}

2. 根據座標計算兩點距離

這個太簡單，跳過

3. 根據弦長求弧長

這個能夠畫個圖，幫助理解：

如今已知弦長 $c$, 半徑 $R$, 要求弧 $r$ 的長度
這很簡單, 只需先求出 $∠\alpha$ 的大小 ：

$$ \alpha = \arcsin(c/2/R)\\ r = 2\alpha \cdot R $$

代碼

const {asin} = Math
const R = 6371

r = asin(c/2/R)*2*R

0b10 最終代碼

/**
 * 獲取兩經緯度之間的距離
 * @param {number} e1 點1的東經, 單位:角度, 若是是西經則爲負
 * @param {number} n1 點1的北緯, 單位:角度, 若是是南緯則爲負
 * @param {number} e2
 * @param {number} n2
 */
function getDistance(e1, n1, e2, n2){
    const R = 6371
    const { sin, cos, asin, PI, hypot } = Math
    
    /** 根據經緯度獲取點的座標 */
    let getPoint = (e, n) => {
        e *= PI/180
        n *= PI/180
        //這裏 R* 被去掉, 至關於先求單位圓上兩點的距, 最後會再將這個距離放大 R 倍
        return {x: cos(n)*cos(e), y: cos(n)*sin(e), z: sin(n)}
    }
    
    let a = getPoint(e1, n1)
    let b = getPoint(e2, n2)
    let c = hypot(a.x - b.x, a.y - b.y, a.z - b.z)
    let r = asin(c/2)*2*R
    return r
}