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

原文链接

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

1、距离

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

hasrsine:

a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ 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位浮点数。据我估计，在简单 余弦公式球面定律 (cosc= cosacosb+ sinasinbcosc)的精度下给出了良好的结果，距离小到地球表面几米。(注意，余弦定律的大地测量形式是规范重新排列的，这样就可以直接使用纬度，而不是用 余纬 )。

这使得简单余弦定理可以在许多大地测量目的（如果不是天文学）中作为 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
       θ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">
    <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>
    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>复制代码