由GPS数据计算地理空间距离

参考链接

https://en.wikipedia.org/wiki/Haversine_formula
https://en.wikipedia.org/wiki/Great-circle_distance
https://blog.csdn.net/u011001084/article/details/52980834
https://www.cnblogs.com/softfair/p/lat_lon_distance_bearing_new_lat_lon.html
https://blog.csdn.net/weixin_33895604/article/details/93930991
https://www.cnblogs.com/osnosn/p/14505778.html
https://blog.csdn.net/allon19/article/details/41649047/
https://blog.csdn.net/mygisforum/article/details/13295223

GPS数据格式

纬度参数格式为：ddmm.mmmmmmm
经度参数格式为：dddmm.mmmmmmm
在实际计算中需先把参数统一为度单位dd.dddddddddd(度)，才能进行后续的转换弧度等操作。
转换算法:
纬度 = dd.0 + (mm.mmmm/60)
经度 = ddd.0 + (mm.mmmm/60)
Int(f / 100) + Frac(f / 100) * 100 / 60;//Frac函数返回实数的小数部分

Haversine方法

用$\lambda ,\varphi ,\Delta \lambda ,\Delta \varphi$表示两个点的地理经度和纬度（弧度）以及他们的绝对差异，$\Delta \sigma$为它们之间的中心角，由余弦的球法则给出。

找到中心角$\Delta \sigma$。给定该角度以弧度为单位，半径为r的球体上的实际弧长d可以用以下公式计算：

Great-circle distance球面余弦公式中有cos(jb-ja)项，当系统的浮点运算精度不高时，在求算较近的两点间的距离时会有较大误差，如果两点在地球表面相距一公里，则中心角的余弦值接近0.99999999会导致较大的舍入误差。但是对于现代的64位浮点数，上面给出的余弦公式的球形定律对于地球表面大于几米的距离没有严重的舍入误差。Haversine方法消除了cos(jb-ja)项，因此不存在短距离求算时对系统计算精度的过多顾虑的问题，采用了正弦函数，即使距离很小，也能保持足够的有效数字。

• R为地球半径，可取平均值 6371km；

//Haversine方法
function GetDistance(lat1, lng1, lat2, lng2) 
{
var radLat1 = lat1 * Math.PI / 180.0;//将角度转换为弧度
var radLat2 = lat2 * Math.PI / 180.0;
var a = radLat1 - radLat2;
var b = lng1 * Math.PI / 180.0 - lng2 * Math.PI / 180.0;

var s = 2 * Math.asin(Math.sqrt(Math.pow(Math.sin(a / 2), 2) +Math.cos(radLat1) * Math.cos(radLat2) * Math.pow(Math.sin(b / 2), 2)));

s = s * 6378.137; // EARTH_RADIUS;
s = Math.round(s * 10000) / 10000;//取整：+0.5向下取整
return s;
            }

//Haversine方法
public static double distHaversineRAD(double lat1, double lon1, double lat2, double lon2) 
{
         double hsinX = Math.sin((lon1 - lon2) * 0.5);
         double hsinY = Math.sin((lat1 - lat2) * 0.5);
         double h = hsinY * hsinY + (Math.cos(lat1) * Math.cos(lat2) * hsinX * hsinX);
         return  2 * Math.atan2(Math.sqrt(h), Math.sqrt(1 - h)) * 6367000;
}

atan2 方法返回一个 -pi 到 pi 之间的数值，表示点 (x, y) 对应的偏移角度。这是一个逆时针角度，以弧度为单位，正X轴和点 (x, y) 与原点连线 之间

因为atan2返回的是弧度值，也就是从-PI到PI,如下图所示，一个半圆是180度=弧度PI,所以1度 = PI/180。

比如现在某个点的坐标为{x:5,y:5},用atan2计算出来的角度degree= Math.atan2(5,5) / (Math.PI/180) 等于45°，注意:这里的第一个参数是y的坐标
但是现在这个角度我们还不能直接使用，因为弧度是一个逆时针方向计算出来的，而我们旋转的时候是按正时针方向旋转，所以我们用的时候要先进行取反:degree = -degree

asin方法返回x 的反正弦值。返回的值是 -PI/2 到 PI/2 之间（即-180度 到 180度）的弧度值。

墨托卡投影(Mercator projection)

（目前仍然不知道咋用）

​ 墨卡托投影是一种“等角正切圆柱投影”，荷兰地图学家墨卡托（Mercator）在1569年拟定：假设地球被围在一个中空的圆柱里，其赤道与圆柱相接触，然后再假想地球中心有一盏灯，把球面上的图形投影到圆柱体上，再把圆柱体展开，这就是一幅标准纬线为零度（即赤道）的“墨卡托投影”绘制出的世界地图。

墨卡托投影在今天对于航海事业起着极为重要的作用，目前世界各国绘制海洋地图时仍广泛使用墨卡托投影，国际水路局（IHB）规定：“除特殊情况外，各国都要用墨卡托投影绘制海图”。国际水路局发行的《大洋水深总图》是把全世界分成24幅编辑的，在南北纬72度之间就是使用墨卡托投影绘成的。

墨卡托投影以整个世界范围，赤道作为标准纬线，本初子午线作为中央经线，两者交点为坐标原点，向东向北为正，向西向南为负。南北极在地图的正下、上方，而东西方向处于地图的正右、左。由于墨卡托投影在两极附近是趋于无限值，因此它并没完整展现了整个世界，地图上最高纬度是85.05度（通过纬度取值范围ys反解计算可得到纬度值为85.05112877980659）。

bool mercator_proj(double B0, double L0, double B, double L, double &X, double&Y)
{
    //标准维度B0，原点维度0，原点经度L0
    static double _A = 6378137, _B = 6356752.3142, _B0 = B0, _L0 = L0;//_B0 = 22 * M_PI / 180, _L0 = 0;
    static double e = sqrt(1 - (_B/_A)*(_B/_A));
    static double e_ = sqrt((_A/_B)*(_A/_B) - 1);
    static double NB0 = ((_A*_A)/_B) / sqrt(1+e_*e_*cos(_B0)*cos(_B0));
    static double K = NB0 * cos(_B0);
    static double E = exp(1);
    
    if(L<-M_PI || L>M_PI || B<-M_PI_2 || B>M_PI_2)
    {
        return false;
    }
    //墨卡托投影正解公式(由B,L得到X,Y——对应的投影坐标)
    Y = K * (L - _L0);
    X = K * log(tan(M_PI_4+B/2) * pow((1-e*sin(B))/(1+e*sin(B)), e/2));
    
    return true;
}