Reference link

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/

GPS data format

The format of latitude parameter is: ddmm mmmmmmm
The format of longitude parameter is: dddmm mmmmmmm
In the actual calculation, the parameters need to be unified into the degree unit dd.dddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd.
Conversion algorithm:
Latitude = dd.0 + (mm.mmmm/60)
Longitude = DDD 0 + (mm.mmmm/60)
Int(f / 100) + Frac(f / 100) * 100 / 60;// The frac function returns the decimal part of a real number

Haversine method

use λ , ϕ , Δ λ , Δ ϕ \lambda , \phi,\Delta\lambda, \Delta\phi λ,ϕ,Δλ,Δϕ Represents the geographic longitude and latitude (radians) of two points and their absolute differences, Δ σ \Delta\sigma Δσ Is the central angle between them, given by the ball law of cosine.

Find the center angle Δ σ \Delta\sigma Δσ. Given that the angle is in radians, the actual arc length d on the sphere with radius r can be calculated by the following formula:

There is cos (JB Ja) term in the spherical cosine formula of great circle distance. When the floating-point operation accuracy of the system is not high, there will be a large error in calculating the distance between the two closer points. If the two points are one kilometer away from the earth's surface, the cosine value of the center angle close to 0.99999 will lead to a large rounding error. However, for modern 64 bit floating-point numbers, the spherical law of the cosine formula given above has no serious rounding error for distances greater than a few meters from the earth's surface. Haversine method eliminates the cos (JB Ja) term, so there is no problem of too much concern about the calculation accuracy of the system in short distance calculation. It adopts sinusoidal function, which can maintain enough significant figures even if the distance is very small.

• R is the radius of the earth, with an average value of 6371km;

//Haversine method
function GetDistance(lat1, lng1, lat2, lng2) 
{
var radLat1 = lat1 * Math.PI / 180.0;//Convert angles to radians
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;//Rounding: + 0.5 rounding down
return s;
            }

//Haversine method
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;
}

The atan2 method returns a value from - pi to PI, representing the offset angle corresponding to the point (x, y). This is a counterclockwise angle, in radians, between the positive X axis and the line between the point (x, y) and the origin

Because atan2 returns radian value, that is, from - pi to PI, as shown in the figure below, a semicircle is 180 degrees = radian PI, so 1 degree = PI/180.

For example, now the coordinates of a point are {x:5,y:5}, and the angle calculated by atan2 is degree = math Atan2 (5,5) / (math. Pi / 180) is equal to 45 °. Note: the first parameter here is the coordinate of y
But now we can't use this angle directly, because the radian is calculated in a counterclockwise direction, and when we rotate, we rotate in a clockwise direction, so when we use it, we need to take the opposite: degree = -degree

The asin method returns the arcsine value of x. The returned value is the radian value between - PI/2 and PI/2 (i.e. - 180 degrees to 180 degrees).

Mercator projection

(I still don't know how to use it)

Mercator projection is a kind of "equiangular tangent cylindrical projection", which was formulated by the Dutch cartographer Mercator in 1569: suppose that the earth is surrounded by a hollow cylinder, and its equator contacts the cylinder, then imagine that there is a lamp in the center of the earth, project the figure on the sphere onto the cylinder, and then expand the cylinder, This is a map of the world drawn by the "Mercator projection" with a standard latitude of zero degrees (i.e. the equator).

Mercator projection plays an extremely important role in navigation today. At present, Mercator projection is still widely used in drawing ocean maps all over the world. The International Waterway Bureau (IHB) stipulates that "except under special circumstances, all countries should use Mercator projection to draw nautical charts". The general map of ocean water depth issued by the International Waterway Bureau is edited by dividing the world into 24 pieces. It is drawn by Mercator projection between North and South latitudes of 72 degrees.

Mercator projection takes the whole world, the equator as the standard latitude and the primary meridian as the central longitude. The intersection of the two is the coordinate origin, which is positive from east to North and negative from west to south. The north and south poles are directly below and above the map, while the east-west direction is directly right and left of the map. Since the Mercator projection tends to be infinite near the poles, it does not fully show the whole world. The highest latitude on the map is 85.05 degrees (the latitude value is 85.05112877980659 can be obtained through the inverse solution of the latitude value range ys).

bool mercator_proj(double B0, double L0, double B, double L, double &X, double&Y)
{
    //Standard dimension B0, origin dimension 0, origin longitude 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;
    }
    //Mercator projection forward solution formula (X,Y -- corresponding projection coordinates are obtained from B and L)
    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;
}