How to calculate the distance between two points on the earth (with formula derivation)

前段时间,Read some algorithms for electronic fences,I'm a little confused about one of the codes that calculates the distance between two points on a sphere,Then I found a related algorithm,The relevant calculation formula is recorded in the great circle distance entry of Wikipedia,The general idea is to find the cosine or sine of the central angle corresponding to the arc length between these two points,Then use the inverse trigonometric function to calculate the radians of the central angle,最后求出：弧长=弧度值 × 地球半径.

注：The above picture uses the Baidu map ranging function,Measure the exit of the railway station in Xiangyang City, Hubei Province and the subway in Changchun City, Jilin Province1The distance from Changchun Railway Station North Subway Station where the line passes

一、具体实现

Suppose there are two points on the sphereA(λ₁ , φ₁)、B(λ₂, φ₂),λ 和 φ respectively represent their longitudes in the map、纬度,θ 为 AB the corresponding central angle,There are two ways to find the radian corresponding to the arc length of two points on the sphere：

$公式1（球面余弦定理）：\theta =acos(cos{\phi }_{1}cos{\phi }_{2}cos({\lambda }_2-{\lambda }_1)+sin{\phi }_{1}sin{\phi }_2)$

$公式2（半正矢量公式）：\theta =2\sqrt{si{n}^2\frac{{\phi }_2-{\phi }_1}{2}+cos{\phi }_{1}cos{\phi }_{2}si{n}^2\frac{{\lambda }_2-{\lambda }_1}{2}}$

JavaScriptThe calculation of trigonometric functions requires the use of radians,The angle needs to be converted to radians first,具体代码实现如下：

var earchRadius = 6371.009;

function toRadians(degrees) {
    
    return degrees * (Math.PI / 180);
}
// 方法一：Use the spherical cosine law
function getDistance(point1, point2) {
    
    var lat1 = toRadians(point1.lat);
    var lat2 = toRadians(point2.lat);
    var lngDiffer = toRadians(point1.lng - point2.lng);

    var rad = Math.acos(Math.cos(lat1) * Math.cos(lat2) * Math.cos(lngDiffer) + Math.sin(lat1) * Math.sin(lat2));

    return rad * earchRadius;
}
// 方法二：Use the half-positive vector formula
function getDistance2(point1, point2) {
    
    var lat1 = toRadians(point1.lat);
    var lat2 = toRadians(point2.lat);
    var latDiffer = toRadians(point1.lat - point2.lat);
    var lngDiffer = toRadians(point1.lng - point2.lng);

    var rad = 2 * Math.asin(Math.sqrt(Math.pow(Math.sin(latDiffer / 2), 2) + Math.cos(lat1) * Math.cos(lat2) * Math.pow(Math.sin(lngDiffer / 2), 2)));

    return rad * earchRadius;
}

注：The average radius of the Earth is approx6,371.009千米,Later in the formula derivation will beR 表示.

下面我们可以验证一下,Let the exit of Xiangyang Railway Station in Xiangyang City, Hubei Province on the map beA点(112.162432,32.064491),A subway in Changchun, Jilin Province1The Changchun Railway Station North Subway Station, which the line passes through, is B点(125.329974,43.919799),The results by the above calculation are as follows：

This matches the distance measured in our opening diagram.

二、公式推导

1、Add auxiliary lines to construct right triangles

从A、BFrom two points, two rays are drawn perpendicular to the equatorial plane,垂足为A₁、B₁;过 A 向 BB₁做垂线,垂足为H.

利用△AHBSolve with the cosine law,具体过程如下：

已知∠BOB₁为Bpoint latitude value,∠AOA₁为Apoint latitude value;∠AOH、∠A₁OB₁为B点和AThe difference between the longitudes of the points,则有

$O{B}_1=R\times cos\angle BO{B}_1,O{A}_1=R\times cos\angle AO{A}_1$

$B{B}_1=R\times sin\angle BO{B}_1,A{A}_1=R\times sin\angle AO{A}_1$

$BH=B{B}_1-H{B}_1=B{B}_1-A{A}_1$

$BH=Rsin\angle BO{B}_1-Rsin\angle AO{A}_1$

$A{H}^2=A_1{B}_1^2$
$A{H}^2=O{B}_1^2+O{A}_1^2-2O{B}_{1}O{A}_{1}cos\angle A_{1}O{B}_1$
$A{H}^2=R^{2}co{s}^2\angle BO{B}_1+R^{2}co{s}^2\angle AO{A}_1-2R^{2}cos\angle BO{B}_{1}cos\angle AO{A}_{1}cos{\lambda }_2-{\lambda }_1$

$A{B}^2=A{H}^2+B{H}^2$

$A{B}^2=A{H}^2+R^{2}si{n}^2\angle BO{B}_1+R^{2}si{n}^2\angle AO{A}_1-2R^{2}sin\angle BO{B}_{1}sin\angle AO{A}_1$

$A{B}^2=2R^2-2R^2(sin\angle BO{B}_{1}sin\angle AO{A}_1+cos\angle BO{B}_{1}cos\angle AO{A}_{1}cos({\lambda }_2-{\lambda }_1))$

在△AOB中,AB² = R² + R² - 2R²cosθ,故可以计算出 θ 的余弦值：

$cos\theta =sin\angle BO{B}_{1}sin\angle AO{A}_1+cos\angle BO{B}_{1}cos\angle AO{A}_{1}cos△\lambda$

即

$cos\theta =cos{\phi }_{1}cos{\phi }_{2}cos({\lambda }_2-{\lambda }_1)+sin{\phi }_{1}sin{\phi }_2$

This is also the cosine law of spherical triangles.

---

2、The distance between two points in a 3D coordinate system

Create a coordinate system as shown in the figure,其中ZThe axis coincides with the earth axis,XOYThe plane is the plane where the equator is located,XOZThe plane is the plane of the prime meridian,则可计算：

$A为（Rcos{\phi }_{1}cos{\lambda }_1,Rcos{\phi }_{1}sin{\lambda }_1,Rsin{\phi }_1）$

$B为（Rcos{\phi }_{2}cos{\lambda }_2,Rcos{\phi }_{2}sin{\lambda }_2,Rsin{\phi }_2）$

$△X=X_B-X_A=Rcos{\phi }_{2}cos{\lambda }_2-Rcos{\phi }_{1}cos{\lambda }_1$

$△Y=Y_B-Y_A=Rcos{\phi }_{2}sin{\lambda }_2-Rcos{\phi }_{1}sin{\lambda }_1$

$△Z=Z_B-Z_A=Rsin{\phi }_2-Rsin{\phi }_1$

$|AB|=\sqrt{△Z^2+△Y^2+△X^2}$

故 $cos\theta =\frac{|OA{|}^2+|OB{|}^2-|AB{|}^2}{2|OA||OB|}$ ,即 $cos\theta =1-\frac{△Z^2+△Y^2+△X^2}{2R^2}$

化简后的结果为

$cos\theta =cos{\phi }_{1}cos{\phi }_{2}cos({\lambda }_2-{\lambda }_1)+sin{\phi }_{1}sin{\phi }_2$

注：This method is briefly described in Wikipedia,The final simplification result is also the spherical cosine theorem

---

3、The haversine formula for a spherical triangle

Just getting started with the semi-positive vector formula,我也是一脸懵.There is no derivation about this formula in a lot of online searches,Until I saw an article on Douban半正矢公式,The semi-positive vector formula is deduced in the article using geometrical methods,Then someone in the comment area pointed out that the haversine function can be substituted into the spherical cosine theorem,我这才恍然大悟.In fact, this formula is not difficult,Still using the half-width formula at all.

Here we need to use the sine vector that is basically not used in trigonometric functions、Haversine function：

$versin\theta =2si{n}^2(\frac{\theta }{2})=1-cos\theta (正矢)$

$haversin\theta =\frac{versin\theta }{2}=\frac{1-cos\theta }{2}(半正矢)$

From this, the equation can be obtained：$si{n}^2(\frac{\theta }{2})=\frac{1-cos\theta }{2}$

然后,Calculated using the previous two methods cosθ Substitute the value into the equation：

$si{n}^2(\frac{\theta }{2})=\frac{1-1+\frac{△Z^2+△Y^2+△X^2}{2R^2}}{2}$

由于 R² Can be removed after expansion,The following calculations are not written R² 了.

$si{n}^2(\frac{\theta }{2})=\frac{△Z^2+△Y^2+△X^2}{4}$

$si{n}^2(\frac{\theta }{2})=\frac{2-2sin{\phi }_{1}sin{\phi }_2-2cos{\phi }_{1}cos{\phi }_{2}cos({\lambda }_2-{\lambda }_1)}{4}（公式1）$

Here we use the product and difference formula in trigonometric functions：

$sin{\phi }_{1}sin{\phi }_2=\frac{cos({\phi }_2-{\phi }_1)-cos({\phi }_2+{\phi }_1)}{2}$

Continue to expand using the square and formula $cos({\phi }_2+{\phi }_1)$ 得：

$sin{\phi }_{1}sin{\phi }_2=\frac{cos({\phi }_2-{\phi }_1)+sin{\phi }_{2}sin{\phi }_1-cos{\phi }_{2}cos{\phi }_1)}{2}$

代入公式1中,可得：

$si{n}^2(\frac{\theta }{2})=\frac{2-cos({\phi }_2-{\phi }_1)+cos{\phi }_{1}cos{\phi }_2-sin{\phi }_{1}sin{\phi }_2-2cos{\phi }_{1}cos{\phi }_{2}cos({\lambda }_2-{\lambda }_1)}{4}$

$si{n}^2(\frac{\theta }{2})=\frac{2-cos({\phi }_2-{\phi }_1)-cos{\phi }_{1}cos{\phi }_2-sin{\phi }_{1}sin{\phi }_2+2cos{\phi }_{1}cos{\phi }_2[1-cos({\lambda }_2-{\lambda }_1)]}{4}$

$si{n}^2(\frac{\theta }{2})=\frac{2-cos({\phi }_2-{\phi }_1)-(cos{\phi }_{1}cos{\phi }_2+sin{\phi }_{1}sin{\phi }_2)+2cos{\phi }_{1}cos{\phi }_2[1-cos({\lambda }_2-{\lambda }_1)]}{4}$

$si{n}^2(\frac{\theta }{2})=\frac{2-2cos({\phi }_2-{\phi }_1)+2cos{\phi }_{1}cos{\phi }_2[1-cos({\lambda }_2-{\lambda }_1)]}{4}$

$si{n}^2(\frac{\theta }{2})=\frac{1-cos({\phi }_1-{\phi }_2)}{2}+cos{\phi }_{1}cos{\phi }_2\frac{1-cos({\lambda }_2-{\lambda }_1)}{2}$

Finally, use the half-width formula：

$si{n}^2(\frac{\theta }{2})=si{n}^2\frac{{\phi }_2-{\phi }_1}{2}+cos{\phi }_{1}cos{\phi }_{2}si{n}^2\frac{{\lambda }_2-{\lambda }_1}{2}$

三、总结一下

总的来说,Calculating the distance between two points on the sphere is probably the two formulas above：Spherical cosine law and half-positive vector formula（That is the half angle solution）.Using the certificate vector formula can reduce the error of spherical cosine theorem calculation,That is, in order to eliminate the error of the mathematical function when calculating the cosine.The included angle corresponding to the arc length is 0时,The value in radians is1,The result obtained by the arc cosine function is not0,but a close one0的浮点数.