根据2个经纬度点，计算这2个经纬度点之间的距离(通过经度纬度得到距离)

转载: https://www.cnblogs.com/softfair/p/distance_of_two_latitude_and_longitude_points.html

球面上任意两点之间的距离计算公式可以参考维基百科上的下述文章。

值得一提的是，维基百科推荐使用Haversine公式，理由是Great-circle distance公式用到了大量余弦函数，而两点间距离很短时（比如地球表面上相距几百米的两点），余弦函数会得出0.999...的结果，会导致较大的舍入误差。而Haversine公式采用了正弦函数，即使距离很小，也能保持足够的有效数字。以前采用三角函数表计算时的确会有这个问题，但经过实际验证，采用计算机来计算时，两个公式的区别不大。稳妥起见，这里还是采用Haversine公式。

其中

R为地球半径，可取平均值 6371km；
φ1, φ2 表示两点的纬度；
Δλ 表示两点经度的差值。

代码

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define PI 3.1415926
#define EARTH_RADIUS 6371393

double deg_to_rad(double deg) {
    return deg / 180 * PI;
}

double hav(double theta) {
    double s = sin(theta / 2);
    return s * s;
}

double get_distance(double lat1, double lng1, double lat2, double lng2) {
    lat1 = deg_to_rad(lat1);
    lng1 = deg_to_rad(lng1);
    lat2 = deg_to_rad(lat2);
    lng2 = deg_to_rad(lng2);

    double dlng = fabs(lng1 - lng2);
    double dlat = fabs(lat1 - lat2);

    double h = hav(dlat) + cos(lat1) * cos(lat2) * hav(dlng);
    double distance = 2 * EARTH_RADIUS * asin(sqrt(h));

    return distance;
}

公式推导

VERSINE(F) = 1 - cos(F)

Haversine名字来历是Ha-VERSINE，即Half-Versine ，表示sin的一半的意思。

hav(A) = (1 - cos(A))/2 = sin(A/2) * sin(A/2)

推倒过程：

如下一个半径为1的圆，O是圆心，A、B是弦(chord)。角度AOB=theta。则角度AOC=theta/2。OC是垂直于AB的垂线(perpendicular)。

AC = sin(theta/2)
AB = 2 * sin(theta/2)

如下地球图所示，假设半径R为1，O是球心，A (lat1,lon1) 和 B (lat2,lon2) 是我们感兴趣的2个点。2跟经度线 lon1，lon2相交于北极(north pole)N。EF所在的线是赤道(equator)。ACBD是平面上的等腰梯形的四个顶点(vertice)。AC和DB的弦（直线）在图上没有画出。CD的位置是：C (lat2,lon1) and D (lat1,lon2)。角度AOC是A点与C点的纬度差 dlat。角度EOF是经度E点和经度F点的差dlon。

弦AC的长度，参照图1的方式，那么是

AC = 2 * sin(dlat/2)

BD也是一样的长度

BD = 2 * sin(dlat/2)

E, F 2个点是赤道上的2个点，它们的纬度是0。EF的距离是

EF = 2 * sin(dlon/2)

A, D 2个点所在的纬度是lat1。AD所在纬度的圆平面的半径是cos(lat1)。从A作一条垂线(perpendicular)到OE为AG，AO是球半径，则OG=cos(lat1)，即A、D所在纬度圆圈的半径（AO`）。

这时候，AD的弦长

AD = 2 * sin(dlon/2) * cos(lat1)

类似的可以推出CB的长度

CB = 2 * sin(dlon/2) * cos(lat2)

下面看一下如何求AB的长度，回到平面等腰梯形，如下图：

AH是到CB的垂线(perpendicular)，CH= (CB-AD)/2。

根据勾股定理(Pythagorean theorem): 【^2表示2的平方】

AH^2 = AC\(^2\) - CH\(^2\)
= AC\(^2\) - (CB - AD)\(^2\)/4

HB 的长度是HB=AD+CH = AD+(CB-AD)/2 = (CB+AD)/2，根据勾股定理得到：

AB^2 = AH\(^2\) + HB\(^2\)
= AC\(^2\) - (CB - AD)\(^2\)/4 + (CB + AD)\(^2\)/4
= AC\(^2\) + CB * AD

根据前面球面上的求经纬距离的方式，我们已经得到 AC、AD和CB的长度，代入公式得到：

AB^2 = 4 * sin\(^2\)(dlat/2) + 4 * cos(lat1) * cos(lat2) * sin\(^2\)(dlon/2)

假设中间值h 是AB长度一半的平方，如下

h = (AB/2)\(^2\)
= sin\(^2\)(dlat/2) + cos(lat1) * cos(lat2) * sin\(^2\)(dlon/2)

（请参看代码里的h）

最后一步，是求得代表AB长度的弧度AOB。参照图1的方式，我们可以知道

弧度AOB = 2 * arcsin(AC)
= 2 * arcsin(AB/2)
= 2 * arcsin(\(\sqrt{h}\))

计算弧长AB

弧长AB = EARTH_RADIUS * 弧度AOB
= 2 * EARTH_RADIUS * arcsin(\(\sqrt{h}\))

根据2个经纬度点，计算这2个经纬度点之间的距离(通过经度纬度得到距离)

代码

公式推导

推荐阅读