(譯)計算距離、方位和更多經緯度之間的點

計算距離、方位和更多經緯度之間的點。最近在研究預測將來座標和速度、時間之間的關係,但願這篇文章對地圖應用有所幫助。

原文連接javascript


該頁面介紹了對經緯度點的各類計算,以及實現它們的公式和代碼片斷。

全部這些公式都是基於球形地球的計算(忽略橢球效應) - 大多數狀況下都是準確的(事實上,地球是很是略微橢圓形的, 使用球形模型的偏差一般高達0.3%1 - 有關詳細信息,請參閱註釋)。php

一、距離

這使用' hasrsine '公式1來計算兩點之間的大圓弧距離 - 即地球表面上的最短距離 - 給出點之間的最短距離。html

hasrsine: java

 a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2) git

 公式:github

 c = 2 ⋅ atan2( √a, √(1−a) )d = R ⋅ cweb

其中φ爲緯度,λ爲經度,R爲地球半徑(平均半徑= 6,371km)。注意角度必須是弧度才能傳遞給三角函數!編程

JavaScript:    
var R = 6371e3; // metres
var φ1 = lat1.toRadians();
var φ2 = lat2.toRadians();
var Δφ = (lat2-lat1).toRadians();
var Δλ = (lon2-lon1).toRadians();

var a = Math.sin(Δφ/2) * Math.sin(Δφ/2) +
        Math.cos(φ1) * Math.cos(φ2) *
        Math.sin(Δλ/2) * Math.sin(Δλ/2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));

var d = R * c;複製代碼

在這些腳本中,我一般用 lat/lon 表示度,以φ/λ表示的距離/長度。發現混合度和半徑是常常讓人抓耳撓腮的問題。bash

即便在很小的距離條件下,「haversine」公式也能夠很好的進行不一樣於基於餘弦球面定律的數值計算。app

以上使用的"(從新)正矢"是 1−cosθ ,以及 「半正矢」是 (1−cosθ)/2 或者 sin²(θ/2)。曾被航海家普遍使用,羅傑·辛諾特於1984年在《天空與望遠鏡》雜誌上對其進行了描述:辛諾特解釋說,大熊星座中Mizar(大熊星座中的一個恆星)和Alcor(北斗七星中的第五等星)之間的角間距(0°11’49.69」)能夠經過使用haversine在TRS-80上精確計算出來。

奇怪的是,c是弧度的角距離,a是兩點弦長一半的平方。

若是 atan2 不可用, c 能夠從 2 ⋅ asin( min(1, √a) ) 計算得出(包括防止四捨五入偏差)。

在中等核的i5我的電腦上使用Chrome,距離計算大約須要2 - 5微秒(所以大約每秒200,000 - 500,000次)。經過公因子分解幾乎沒有用,由於JIT編譯器可能會對它們進行優化。

二、餘弦球面定律

事實上,JavaScript(以及大多數現代計算機和語言)使用提供了15個重要的精度數字的 「IEEE 754」 64位浮點數。據我估計,在簡單餘弦公式球面定律(cosc= cosacosb+ sinasinbcosc)的精度下給出了良好的結果,距離小到地球表面幾米。(注意,餘弦定律的大地測量形式是規範從新排列的,這樣就能夠直接使用緯度,而不是用餘緯)。

這使得簡單餘弦定理能夠在許多大地測量目的(若是不是天文學)中做爲 haversine 公式的另外一個選擇。這種選擇多是由編程語言、處理器、編碼上下文、可用地三角函數(在不一樣語言中)等驅動的,而且,在很是小地距離中。等矩形的近似值可能更合適。

餘弦定理 :d = acos( sin φ1 ⋅ sin φ2 + cos φ1 ⋅ cos φ2 ⋅ cos Δλ ) ⋅ R

JavaScript:    
var φ1 = lat1.toRadians(), φ2 = lat2.toRadians(), Δλ = (lon2-lon1).toRadians(), R = 6371e3; // gives d in metres
var d = Math.acos( Math.sin(φ1)*Math.sin(φ2) + Math.cos(φ1)*Math.cos(φ2) * Math.cos(Δλ) ) * R;複製代碼
Excel:  
=ACOS( SIN(lat1)*SIN(lat2) + COS(lat1)*COS(lat2)*COS(lon2-lon1) ) * 6371000

(or with lat/lon in degrees): 
=ACOS( SIN(lat1*PI()/180)*SIN(lat2*PI()/180) + COS(lat1*PI()/180)*COS(lat2*PI()/180)*COS(lon2*PI()/180-lon1*PI()/180) ) * 6371000複製代碼

雖然比較簡單,但在個人測試中,餘弦定理比 haversine 略慢。

三、等矩形近似值

若是性能是一個問題,準確性就不那麼重要了。在很小距離上,勾股定理能夠在這個等矩形投影中使用:

公式:

 x = Δλ ⋅ cos φm

y = Δφ

d = R ⋅ √x² + y²

JavaScript:    
var x = (λ21) * Math.cos((φ12)/2);
var y = (φ21);
var d = Math.sqrt(x*x + y*y) * R;複製代碼

這僅僅使用一個三角函數和一個根號函數,而cos定律中有6個三角函數,而 haversine 中有7個三角函數+ 2根號函數。精度有點複雜:沿着子午線沒有偏差,不然它們依賴於距離、方位和緯度,但對於許多用途來講足夠小(一般與球面近似自己相比微不足道)。

此外,還可使用極座標平面地球公式:使用餘緯 θ1 = π/2−φ1 and θ2 = π/2−φ2,而後 θ1 = π/2−φ1 and θ2 = π/2−φ2。我沒有 比較精度。

四、方位

通常來講,當你沿着一個大圓弧路徑(直角方向)移動時,當前朝向會發生變化。根據距離和緯度的不一樣,最終的航向與最初的航向會有不一樣程度的差別(好比從35°N 45°E(≈Baghdad)到35°N 135°E(≈Osaka),你會從一個60°的航向開始,最終到達一個120°的航向)。

這個公式適用於若是初始方位(有時稱爲前方位角)沿着大圓弧走一條直線,你會從起點走到終點1

公式: θ = atan2( sin Δλ ⋅ cos φ2 , cos φ1 ⋅ sin φ2 − sin φ1 ⋅ cos φ2 ⋅ cos Δλ )

φ11是起點,φ22是終點(Δλ是經度的差別)

JavaScript:
(all angles in radians)
var y = Math.sin(λ2-λ1) * Math.cos(φ2);
var x = Math.cos(φ1)*Math.sin(φ2) -
        Math.sin(φ1)*Math.cos(φ2)*Math.cos(λ2-λ1);
var brng = Math.atan2(y, x).toDegrees();複製代碼
Excel:
(all angles in radians)
=ATAN2(COS(lat1)*SIN(lat2)-SIN(lat1)*COS(lat2)*COS(lon2-lon1), 
       SIN(lon2-lon1)*COS(lat2)) 複製代碼

注意,Excel將參數轉換爲ATAN2—請參閱下面的註釋

因爲 atan2 返回值範圍中的值-π……+π(-180°…+180°),將結果正常化爲羅盤方位(範圍爲0°…360°,- ve值轉換爲180°範圍…360°),轉換爲度,而後使用(θ+ 360)% 360,%是(浮點數)模。

對最後的方位,只是把初始方位從終點到起點並反轉它(使用θ=θ+ 180)% 360)。

五、中點

這是兩點之間沿大圓路徑的中點1

公式: 

 Bx = cos φ2 ⋅ cos Δλ

By = cos φ2 ⋅ sin Δλ

φm = atan2( sin φ1 + sin φ2, √(cos φ1 + Bx)² + By² )
λm = λ1 + atan2(By, cos(φ1)+Bx)

JavaScript:
(all angles in radians)
var Bx = Math.cos(φ2) * Math.cos(λ2-λ1);
var By = Math.cos(φ2) * Math.sin(λ2-λ1);
var φ3 = Math.atan2(Math.sin(φ1) + Math.sin(φ2),
              Math.sqrt( (Math.cos(φ1)+Bx)*(Math.cos(φ1)+Bx) +By*By ) );
var λ3 = λ1 + Math.atan2(By, Math.cos(φ1) + Bx);複製代碼

經度可使用(lon+540)%360-180歸一化爲-180…+180

正如初始方位可能與最終方位不一樣,中點也可能不在經緯度之間。35°N、45°E和35°N、135°E之間的中點在45°N、90°E附近。

六、中間點

也能夠計算出兩點之間沿大圓弧路徑上任意處的中間點1

公式:

a = sin((1−f)⋅δ) / sin δ

b = sin(f⋅δ) / sin δ

x = a ⋅ cos φ1 ⋅ cos λ1 + b ⋅ cos φ2 ⋅ cos λ2

y = a ⋅ cos φ1 ⋅ sin λ1 + b ⋅ cos φ2 ⋅ sin λ2

z = a ⋅ sin φ1 + b ⋅ sin φ2

φi = atan2(z, √x² + y²)

λi = atan2(y, x)

f是沿大圓航線的一部分(f = 0是點1,f = 1是點2),δ是兩點之間的角距離d / R。

七、目標點給出從起點到終點的距離和方位方位

已知一個起始點、初始方位和距離,這將計算沿着(最短距離)大圓弧運動的終點和最終方位。

公式:φ2 = asin( sin φ1 ⋅ cos δ + cos φ1 ⋅ sin δ ⋅ cos θ )

λ2 = λ1 + atan2( sin θ ⋅ sin δ ⋅ cos φ1, cos δ − sin φ1 ⋅ sin φ2 )

φ是緯度,λ是經度,θ是軸承方位(順時針從北),δ是角距離d / R,d是通過的距離,R是地球的半徑

JavaScript:
(all angles in radians)
var φ2 = Math.asin( Math.sin(φ1)*Math.cos(d/R) +
                    Math.cos(φ1)*Math.sin(d/R)*Math.cos(brng) );
var λ2 = λ1 +Math.atan2(Math.sin(brng)*Math.sin(d/R)*Math.cos(φ1),
                         Math.cos(d/R)-Math.sin(φ1)*Math.sin(φ2));複製代碼

經度可使用 (lon+540)%360-180 規範化爲 -180…+180

Excel:
(all angles in radians)
lat2: =ASIN(SIN(lat1)*COS(d/R) + COS(lat1)*SIN(d/R)*COS(brng))
lon2: =lon1 + ATAN2(COS(d/R)-SIN(lat1)*SIN(lat2), SIN(brng)*SIN(d/R)*COS(lat1)) * Remember that Excel reverses the arguments to ATAN2 – see notes below複製代碼

對於最終方位,只需將初始方位從起點移到起點,而後用(brng+180)%360將其反轉。

九、給定起點和方位兩條路徑的交點

這是一個比這一頁上的大多數計算要複雜得多的計算,可是我已經被問過不少次了。這來自埃德·威廉的航空公式。參見下面的JavaScript。

公式: δ12 = 2⋅asin( √(sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)) ) (角距離. p1–p2)

θa = acos( ( sin φ2 − sin φ1 ⋅ cos δ12 ) / ( sin δ12 ⋅ cos φ1 ) ) (初始/ 最終 方位)

θa = acos( ( sin φ2 − sin φ1 ⋅ cos δ12 ) / ( sin δ12 ⋅ cos φ1 ) ) (點1和點2之間)

if sin(λ2−λ1) > 0
       θ12 = θa
       θ21 = 2π − θb
    else
       θ12 = 2π − θa
       θ21 = θb    
       α1 = θ13 − θ12          (角度 p2–p1–p3)
       α2 = θ21 − θ23          (角度 p1–p2–p3)複製代碼

α3 = acos( −cos α1 ⋅ cos α2 + sin α1 ⋅ sin α2 ⋅ cos δ12 ) (角度 p1–p2–p3)

δ13 = atan2( sin δ12 ⋅ sin α1 ⋅ sin α2 , cos α2 + cos α1 ⋅ cos α3 ) (角距離 p1–p3)

φ3 = asin( sin φ1 ⋅ cos δ13 + cos φ1 ⋅ sin δ13 ⋅ cos θ13 ) (p3緯度)

Δλ13 = atan2( sin θ13 ⋅ sin δ13 ⋅ cos φ1 , cos δ13 − sin φ1 ⋅ sin φ3 ) (經度p1-p3)

λ3 = λ1 + Δλ13 (p3 經度)

φ1、λ1、θ13:1起點&(初始)軸承從1日指向交點

φ2、λ2、θ23:2日起點&(初始)軸承2指向交點

φ3、λ3:交點

% =(浮點數)模

若是sin α1= 0 以及 sin α2 = 0:無窮多解,若是sin α1 ⋅ sin α2 < 0:歧義解。這個公式並不老是適用於經向線或赤道線。

用向量比用球面三角簡單多了: latlong-vectors.html.

十、航跡距離

這是一個新的問題:有時我被問到關於一個點到大圓弧路徑的距離(有時稱爲交叉軌跡偏差)。

公式: 

dxt = asin( sin(δ13) ⋅ sin(θ13−θ12) ) ⋅ R

δ13是(角)距離起點第三點

θ13是(初始)軸承從起點到第三點

θ 12是(初始)軸承從起點到終點
R是地球的半徑
JavaScript:    
var δ13 = d13 / R;
var dAt = Math.acos(Math.cos(δ13)/Math.cos(dXt/R)) * R;複製代碼

十一、離兩極最近的點

「克萊勞公式」會用給出的方位 θ 和最大圓上的緯度 φ,計算出一個大圓弧路徑的最大緯度:
公式:φ max = acos( | sin θ ⋅ cos φ | )
JavaScript:    
var φMax = Math.acos(Math.abs(Math.sin(θ)*Math.cos(φ)));複製代碼

十二、恆向線

「恆向線」(或斜航線)是一條恆定方位的路徑,它以相同的角度穿過全部子午線。

水手們過去經常(如今有時仍然)沿着恆向線航行,由於跟隨着一個恆定的羅盤方位比跟隨繞一個大圓須要不斷調整方位要容易得多。等高線是墨卡託投影圖上的直線(也有助於導航)。

等高線通常比大圓線長。例如,從倫敦到紐約,沿一條斜脊線比沿一個大圓要長4%——這對航空燃料很重要,但對帆船來講並不特別重要。從紐約到北京——接近可能的最極端的例子(雖然不能航行)——沿着一條斜脊線要長30%。

計算恆向線的關鍵是範德曼函數,它給出了已知緯度的墨卡託投影地圖上的高度:ln(tanφ + secφ) 或是 ln( tan(π/4+φ/2) )。固然,這在兩極是趨於無窮的(與墨卡託投影一致)。對於強迫症患者,甚至有一個橢圓形的版本,「等量緯度」:ψ= ln (tan(π/ 4 +φ/ 2)/ ((1−e⋅sinφ)/ (1 + e⋅sinφ)]e / 2),或其更好的等值於 ψ= atanh (sinφ)−e⋅atanh (e⋅sinφ)。

該公式從球形緯度和經度座標推算出墨卡託投影以東和以北值,而後:

E = R ⋅ λ

N = R ⋅ ln( tan(π/4 + φ/2) )

如下公式來源於埃德·威廉姆斯的航空公式1

1.距離

因爲恆向線是墨卡託投影上的直線,沿恆向線的兩點之間的距離就是該線的長度(畢達哥拉斯),可是項目的失真須要獲得補償。

按固定緯度航行(東西航行),cosφ 補償很簡單,一般狀況下,它是Δφ/Δψ,Δψ = ln( tan(π/4 + φ2/2) / tan(π/4 + φ1/2) ) (「投影」緯度差)

公式:

Δψ = ln( tan(π/4 + φ2/2) / tan(π/4 + φ1/2) ) (「投影」緯度差)

q = Δφ/Δψ (or cosφ for E-W line)

d = √(Δφ² + q²⋅Δλ²) ⋅ R

φ是緯度,λ是經度,Δλ是正在走的最短路線(< 180°),R是地球的半徑,ln是天然對數

JavaScript:
(all angles in radians)
var Δψ =Math.log(Math.tan(Math.PI/4+φ2/2)/Math.tan(Math.PI/4+φ1/2));
var q = Math.abs(Δψ) > 10e-12 ? Δφ/Δψ : Math.cos(φ1); // E-W course becomes ill-conditioned with 0/0

// if dLon over 180° take shorter rhumb line across the anti-meridian:
if (Math.abs(Δλ) > Math.PI) Δλ = Δλ>0 ? -(2*Math.PI-Δλ) : (2*Math.PI+Δλ);

var dist = Math.sqrt(Δφ*Δφ + q*q*Δλ*Δλ) * R;複製代碼

2.方位

恆向線是墨卡託投影上的一條直線,其投影角度與羅盤方位相等。

公式:

Δψ = ln( tan(π/4 + φ2/2) / tan(π/4 + φ1/2) ) (「投影」緯度差)

θ = atan2(Δλ, Δψ)

φ是緯度,λ是經度,Δλ是正在走的最短路線(< 180°),R是地球的半徑,ln是天然對數

JavaScript:
(all angles in radians)
var Δψ =Math.log(Math.tan(Math.PI/4+φ2/2)/Math.tan(Math.PI/4+φ1/2));

// if dLon over 180° take shorter rhumb line across the anti-meridian:
if (Math.abs(Δλ) > Math.PI) Δλ = Δλ>0 ? -(2*Math.PI-Δλ) : (2*Math.PI+Δλ);

var brng = Math.atan2(Δλ, Δψ).toDegrees();複製代碼

3.終點

給定一個起點和一個距離d和不變的方位θ,這將計算出終點。若是你沿着一條等高線保持恆定的方位,你就會逐漸地向其中一個極點螺旋靠近。

公式:

δ = d/R (角距離)

φ2 = φ1 + δ ⋅ cos θ

Δψ = ln( tan(π/4 + φ2/2) / tan(π/4 + φ1/2) ) (「投影」緯度差)

q = Δφ/Δψ (or cos φ for E-W line)

Δλ = δ ⋅ sin θ / q

λ2 = λ1 + Δλ

φ是緯度,λ是經度,Δλ是正在走的最短路線(< 180°),ln是天然對數,R是地球的半徑

JavaScript:
(all angles in radians)
var δ = d/R;
var Δφ = δ * Math.cos(θ);
var φ2 = φ1 + Δφ;

var Δψ =Math.log(Math.tan(φ2/2+Math.PI/4)/Math.tan(φ1/2+Math.PI/4));
var q = Math.abs(Δψ) > 10e-12 ? Δφ / Δψ : Math.cos(φ1); // E-W course becomes ill-conditioned with 0/0

var Δλ = δ*Math.sin(θ)/q;
var λ2 = λ1 + Δλ;

// check for some daft bugger going past the pole, normalise latitude if so
if (Math.abs(φ2) > Math.PI/2) φ2 = φ2>0 ? Math.PI-φ2 : -Math.PI-φ2;複製代碼

4.中間點

這個計算「斜向中點」的公式是由羅伯特·希爾(Robert Hill)和克萊夫·圖特(Clive Tooth1) (thx Axel!)推導出來的。

公式:

φm = (φ1+φ2) / 2

f1 = tan(π/4 + φ1/2)

f2 = tan(π/4 + φ2/2)

fm = tan(π/4+φm/2)

λm = [ (λ2−λ1) ⋅ ln(fm) + λ1 ⋅ ln(f2) − λ2 ⋅ ln(f1) ] / ln(f2/f1)

φ是緯度,經度λ,ln是天然對數

JavaScript:
(all angles in radians)
if (Math.abs(λ2-λ1) > Math.PI) λ1 += 2*Math.PI; // crossing anti-meridian

var φ3 = (φ1+φ2)/2;
var f1 = Math.tan(Math.PI/4 + φ1/2);
var f2 = Math.tan(Math.PI/4 + φ2/2);
var f3 = Math.tan(Math.PI/4 + φ3/2);
var λ3 = ( (λ2-λ1)*Math.log(f3) + λ1*Math.log(f2) - λ2*Math.log(f1) ) / Math.log(f2/f1);

if (!isFinite(λ3)) λ3 = (λ1+λ2)/2; // parallel of latitude複製代碼

經度可使用(lon+540)%360-180規範化爲-180…+180。

5.在web頁面中使用腳本

在web頁面中使用這些腳本大體以下:

<!doctype html>
<html lang="en">
<head>
    <title>Using the scripts in web pages</title>
    <meta charset="utf-8">
    <script type="module">
        import LatLon from 'https://cdn.jsdelivr.net/gh/chrisveness/geodesy@2.0.0/latlon-spherical.min.js';
        document.addEventListener('DOMContentLoaded', function() {
            document.querySelector('#calc-dist').onclick = function() {
                calculateDistance();
            }
        });
        function calculateDistance() {
            const p1 = LatLon.parse(document.querySelector('#point1').value);
            const p2 = LatLon.parse(document.querySelector('#point2').value);
            const dist = parseFloat(p1.distanceTo(p2).toPrecision(4));
            document.querySelector('#result-distance').textContent = dist + ' metres';
        }
    </script>
</head>
<body>
<form>
    Point 1: <input type="text" name="point1" id="point1" placeholder="lat1,lon1">
    Point 2: <input type="text" name="point2" id="point2" placeholder="lat2,lon2">
    <button type="button" id="calc-dist">Calculate distance</button>
    <output id="result-distance"></output>
</form>
</body>
</html>複製代碼

請參閱下面JavaScript源代碼,也能夠在GitHub上得到。完整的文檔以及測試套件都是可用的。

6.筆記

(1)

精度:因爲地球並不徹底是一個球體,因此在使用球面幾什麼時候存在小的偏差。地球實際上大體是橢球體(或者更精確地說,是扁圓球體),半徑在6,378公里(赤道)和6,357公里(極地)之間變化,局部曲率半徑從6,336公里(赤道子午線)到6,399公里(極地)之間變化。6371公里是地球平均半徑的廣泛接受值。這意味着,根據緯度和旅行方向的不一樣,假設球面幾何形狀穿越赤道的偏差可能高達0.55%,但通常低於0.3% (whuber在stackexchange上對此進行了詳細的研究)。對於我來講,1km中精度超過3m就足夠了,可是若是你想要更高的精度,你可使用Vincenty公式來計算橢球體上的測地線距離,結果能夠精確到1mm之內。(徹底出於反常,我歷來沒有須要過這麼精確的代碼,我查了一下這個公式,發現JavaScript實現比我預想的要簡單)。

(2)

三角函數以弧度爲參數,因此緯度、經度、和方位(十進制或度/分鐘/秒)須要轉換爲弧度,rad =π.deg / 180。當弧度轉換回度時(deg = 180.rad/π),若是使用帶符號的十進制度,則West爲負。對於方位,其在-π到+π(-180°+ 180°)的範圍須要轉換爲0到+ 2π(0°-360°)。這能夠經過(brng+2.π)%2.π(或brng+360)%360),%(浮點數)是模運算符(請注意,不一樣的語言以不一樣的方式實現模操做)。

(3)

全部方位均相對於正北,0°=N, 90°=E等。若是你用指南針工做,地磁北極從正北極以一種複雜的方式圍繞地球發生變化,這種差別必須由本地地圖上顯示的差別來補償。

(4)

這裏普遍使用的atan2()函數接受兩個參數,atan2(y, x),並計算比值y/x的反正切。它比atan(y/x)更靈活,由於它處理x = 0,而且返回在全部4個象限的值-π+π(反正切函數返回在範圍-π/ 2 +π/ 2之間的值)。

(5)

若是你在電子表格(Microsoft Excel, LibreOffice Calc, Google Sheets, Apple Numbers)中使用任何公式包括atan2,你將須要反轉參數,如Excel等,它們與JavaScript相反,常規順序爲atan(y,x),可是Excel使用atan2(x, y)。在(VBA)宏觀中使用 atan2,你能夠用 WorksheetFunction.Atan2()。

(6)

若是您正在使用谷歌地圖,這些函數中的一些(computeDistanceBetween()、computeHeading()、computeOffset()、插值()等)如今在谷歌地圖 API V3「球面」庫中提供。注意,它們使用的默認地球半徑爲6,378,137米)。

(7)

若是你使用英國地形測量網格參考,我已經實現了一個腳本之間的轉換Lat/Long和OS網格參考

(8)

若是您使用UTM座標或MGRS網格引用,我已經實現了用於在Lat/Long、UTM和MGRS之間轉換的腳本。

(9)

我從美國人口普查局的GIS FAQ中學到了不少,如今已經沒有了,因此我複製了一份。

(10)

感謝埃德·威廉姆斯的航空公式。

(11)

對於英里,用km除以1.609344。

(12)

對於海里,用km除以1.852。

7.計算工具

點擊查看

相關文章
相關標籤/搜索