计算距离、方位和更多经纬度之间的点。最近在研究预测未来坐标和速度、时间之间的关系，希望这篇文章对地图应用有所帮助。

作者：狐狸家的鱼

本文链接：计算距离、方位和更多经纬度之间的点

原文链接：Calculate distance, bearing and more between Latitude/Longitude points

GitHub：sueRimn

该页面介绍了对经纬度点的各种计算，以及实现它们的公式和代码片段。

所有这些公式都是基于球形地球的计算（忽略椭球效应） – 大多数情况下都是准确的（事实上，地球是非常略微椭圆形的， 使用球形模型的误差通常高达0.3%¹ – 有关详细信息，请参阅注释）。

1、距离

这使用’ hasrsine ‘公式¹来计算两点之间的大圆弧距离 – 即地球表面上的最短距离 – 给出点之间的最短距离。

                                  hasrsine:     a = sin²(Δφ/2) + cos φ₁ ⋅ cos φ₂ ⋅ sin²(Δλ/2)

            公式:      c = 2 ⋅ atan2( √a, √(1−a) )

d = R ⋅ c

其中φ为纬度，λ为经度，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 表示度，以 φ/λ 表示的距离/长度。发现混合度和半径是经常让人抓耳挠腮的问题。

即使在很小的距离条件下，“haversine”公式也可以很好的进行不同于基于余弦球面定律的数值计算。

以上使用的”(重新)正矢”是 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编译器可能会对它们进行优化。

2、余弦球面定律

事实上，JavaScript(以及大多数现代计算机和语言)使用提供了15个重要的精度数字的 “IEEE 754” 64位浮点数。据我估计，在简单余弦公式球面定律(cos c = cos a cos b + sin a sin b cos C)的精度下给出了良好的结果，距离小到地球表面几米。(注意，余弦定律的大地测量形式是规范重新排列的，这样就可以直接使用纬度，而不是用余纬)。

这使得简单余弦定理可以在许多大地测量目的（如果不是天文学）中作为 haversine 公式的另一个选择。这种选择可能是由编程语言、处理器、编码上下文、可用地三角函数（在不同语言中）等驱动的，并且，在非常小地距离中。等矩形的近似值可能更合适。

余弦定理 ：d = acos( sin φ₁ ⋅ sin φ₂ + cos φ₁ ⋅ cos φ₂ ⋅ 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 略慢。

3、等矩形近似值

如果性能是一个问题,准确性就不那么重要了。在很小距离上，勾股定理可以在这个等矩形投影中使用：

公式：  x = Δλ ⋅ cos φm

y = Δφ

             d = R ⋅ √x² + y²

JavaScript:    
var x = (λ2-λ1) * Math.cos((φ1+φ2)/2);
var y = (φ2-φ1);
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。我没有 比较精度。

4、方位

一般来说，当你沿着一个大圆弧路径(直角方向)移动时，当前朝向会发生变化。根据距离和纬度的不同，最终的航向与最初的航向会有不同程度的差异（比如从35°N 45°E(≈Baghdad)到35°N 135°E(≈Osaka),你会从一个60°的航向开始，最终到达一个120°的航向)。

这个公式适用于如果初始方位(有时称为前方位角)沿着大圆弧走一条直线，你会从起点走到终点^1：

                             公式： θ = atan2( sin Δλ ⋅ cos φ₂ , cos φ₁ ⋅ sin φ₂ − sin φ₁ ⋅ cos φ₂ ⋅ cos Δλ )

φ₁,λ₁是起点,φ₂,λ₂是终点(Δλ是经度的差异)

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)。

5、中点

这是两点之间沿大圆路径的中点¹。

公式：   B_x = cos φ₂ ⋅ cos Δλ

             B_y = cos φ₂ ⋅ sin Δλ

                                                              φ_m = atan2( sin φ₁ + sin φ₂, √(cos φ₁ + B_x)² + B_y² )
                                 λ_m = λ₁ + atan2(B_y, cos(φ₁)+B_x)

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附近。

6、中间点

也可以计算出两点之间沿大圆弧路径上任意处的中间点¹。

公式：a = sin((1−f)⋅δ) / sin δ

    b = sin(f⋅δ) / sin δ

                                              x = a ⋅ cos φ₁ ⋅ cos λ₁ + b ⋅ cos φ₂ ⋅ cos λ₂

                                            y = a ⋅ cos φ₁ ⋅ sin λ₁ + b ⋅ cos φ₂ ⋅ sin λ₂

                  z = a ⋅ sin φ₁ + b ⋅ sin φ₂

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

λ_i = atan2(y, x)

 f是沿大圆航线的一部分(f = 0是点1,f = 1是点2),δ是两点之间的角距离d / R。

7、目标点给出从起点到终点的距离和方位方位

已知一个起始点、初始方位和距离，这将计算沿着(最短距离)大圆弧运动的终点和最终方位。

公式：φ₂ = asin( sin φ₁ ⋅ cos δ + cos φ₁ ⋅ sin δ ⋅ cos θ )

                                λ₂ = λ₁ + atan2( sin θ ⋅ sin δ ⋅ cos φ₁, cos δ − sin φ₁ ⋅ sin φ₂ )

φ是纬度，λ是经度，θ是轴承方位(顺时针从北),δ是角距离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将其反转。

9、给定起点和方位两条路径的交点

这是一个比这一页上的大多数计算要复杂得多的计算，但是我已经被问过很多次了。这来自埃德·威廉的航空公式。参见下面的JavaScript。

公式：     δ₁₂ = 2⋅asin( √(sin²(Δφ/2) + cos φ₁ ⋅ cos φ₂ ⋅ sin²(Δλ/2)) )   （角距离. p1–p2）

                       θ_a = acos( ( sin φ₂ − sin φ₁ ⋅ cos δ₁₂ ) / ( sin δ₁₂ ⋅ cos φ₁ ) )   （初始/ 最终 方位）

                     θ_a = acos( ( sin φ₂ − sin φ₁ ⋅ cos δ₁₂ ) / ( sin δ₁₂ ⋅ cos φ₁ ) )   （点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）

     　　　　　　　　　　　　　　　                     α₃ = acos( −cos α₁ ⋅ cos α₂ + sin α₁ ⋅ sin α₂ ⋅ cos δ₁₂ )   （角度 p1–p2–p3）

              　　　　　　　　　　　　　　　　　                 δ₁₃ = atan2( sin δ₁₂ ⋅ sin α₁ ⋅ sin α₂ , cos α₂ + cos α₁ ⋅ cos α₃ )   （角距离 p1–p3）

  　　　　　　　　　　　　　　　　           φ₃ = asin( sin φ₁ ⋅ cos δ₁₃ + cos φ₁ ⋅ sin δ₁₃ ⋅ cos θ₁₃ )   （p3纬度）

                                                                                  　　　Δλ₁₃ = atan2( sin θ₁3 ⋅ sin δ₁₃ ⋅ cos φ₁ , cos δ₁₃ − sin φ₁ ⋅ sin φ₃ )   （经度p1-p3）

λ₃ = λ₁ + Δλ₁₃   (p3 经度)

                                  φ₁、λ₁、θ₁₃：1起点&(初始)轴承从1日指向交点

                                 φ₂、λ₂、θ₂₃：2日起点&(初始)轴承2指向交点

φ₃、λ₃：交点

% =(浮点数)模

如果sin α₁= 0 以及 sin  α₂ = 0：无穷多解，如果sin α₁ ⋅ sin α₂ < 0：歧义解。这个公式并不总是适用于经向线或赤道线。

用向量比用球面三角简单多了： latlong-vectors.html.

10、航迹距离

这是一个新的问题：有时我被问到关于一个点到大圆弧路径的距离(有时称为交叉轨迹误差)。

　　　　　　 公式： d_xt = asin( sin(δ₁₃) ⋅ sin(θ₁₃−θ₁₂) ) ⋅ R

　　　　δ_13是(角)距离起点第三点

　　　　　　　θ₁3是(初始)轴承从起点到第三点

JavaScript:    
var δ13 = d13 / R;
var dAt = Math.acos(Math.cos(δ13)/Math.cos(dXt/R)) * R;

JavaScript:    
var φMax = Math.acos(Math.abs(Math.sin(θ)*Math.cos(φ)));

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;

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();

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;

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

<!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>