求通过经纬度计算地球表面任意两点间距离的公式

求通过经纬度计算地球表面任意两点间距离的公式
/// <summary>
    /// //经纬度类
    /// vp:hsg
    /// create date:2011-05-04
    /// </summary>
    public class LatLngClass
    {
        public double Lat = 0.0;
        public double Lng = 0.0;
        public LatLngClass()
        {
        }
        public LatLngClass(double pLat,double pLng)
        {
            Lat = pLat;
            Lng = pLng;
        }
        //获取两个经纬度坐标直接的距离
        public double GetDistance(LatLngClass BlatLng)
        {
            double rbc = 0;
            LatLngClass AlatLng = this;
            rbc = GetDistance(AlatLng, BlatLng);
            return rbc;
        }
        //获取两个经纬度坐标直接的距离
        public  double GetDistance(LatLngClass AlatLng,LatLngClass BlatLng)
        {
            double rbc = 0;
            double xgA=0;double ygA=0;
            BlToGs(ref xgA, ref ygA, AlatLng.Lng, AlatLng.Lat);
            double xgB = 0; double ygB = 0;
            BlToGs(ref xgB, ref ygB, BlatLng.Lng, BlatLng.Lat);
            rbc = GetDistance(xgA, ygA, xgB, ygB);
            
            return rbc;
        }
        //经纬度转为高斯平面直角坐标(bl单位度)  
        public void BlToGs(ref double xg,ref double yg, double lon, double lat)         
        {
            int dh;//带号   
            double temp1, temp2, a, b, c, d, e, f;
            double l1, l2;
            int n;
            n = (int)(lon / 6);
            dh = n + 1;
            l2 = 3 + 6 * n;
            l1 = lon - l2;
            temp1 = lat * Math.PI / 180;
            temp2 = l1 * Math.PI / 180;
            a = 6367558.497 * temp1;
            b = (16036.480 - 1597237.956 * Math.Pow(temp2, 2) - 268563.280 * Math.Pow(temp2, 4)) * Math.Sin(2 * temp1);
            c = (16.828 - 1340.831 * Math.Pow(temp2, 2) + 201450.536 * Math.Pow(temp2, 4)) * Math.Sin(4 * temp1);
            xg = a - b + c;
            d = (6383594.975 + 535998.795 * Math.Pow(temp2, 2) + 54206.791 * Math.Pow(temp2, 4)) * Math.Cos(temp1);
            e = (5356.713 - 534204.967 * Math.Pow(temp2, 2) - 134966.691 * Math.Pow(temp2, 4)) * Math.Cos(3 * temp1);
            f = (6.744 + 81276.496 * Math.Pow(temp2, 4)) * Math.Cos(5 * temp1);
            yg = (temp2) * (d - e + f) + 500000 + dh * 1000000;
        }
        //获取直角坐标系统距离
        public  double GetDistance(double Vx,double Vy,double mX,double mY)
        {
            return Math.Sqrt((mX - Vx) * (mX - Vx) + (mY - Vy) * (mY - Vy));
        }
    }
测试方法
static public void LatLng()
        {
            double xg = 0; double yg = 0;
            LatLngClass ll = new LatLngClass(103.021144,30.234332);
            ll.BlToGs(ref xg, ref yg, ll.Lng, ll.Lat);
            Console.WriteLine("x,y=" + xg + "," + yg);
            double x = 0; double y = 0;
            LatLngClass Lc = new LatLngClass(103.021234, 30.234323);
            Lc.BlToGs(ref x, ref y, Lc.Lng, Lc.Lat);
            Console.WriteLine("x,y=" + x + "," + y);
            Console.WriteLine("distance="+ll.GetDistance(Lc));
        }

通过经纬度计算两点间距离的公式

 

对上面的公式解释如下：

1.Lat1 Lung1 表示A点经纬度，Lat2 Lung2 表示B点经纬度；

2.a=Lat1 – Lat2 为两点纬度之差  b=Lung1 -Lung2 为两点经度之差；

3.6378.137为地球半径，单位为公里；

计算出来的结果单位为公里。 

注：公式中经纬度均用弧度表示； 

 

地球是一个近乎标准的椭球体，它的赤道半径为6378.140千米，极半径为6356.755千米，平均半径6371.004千米。如果我们假设地球是一个完美的球体，那么它的半径就是地球的平均半径，记为R。如果以0度经线为基准，那么根据地球表面任意两点的经纬度就可以计算出这两点间的地表距离（这里忽略地球表面地形对计算带来的误差，仅仅是理论上的估算值）。设第一点A的经纬度为(LonA, LatA)，第二点B的经纬度为(LonB, LatB)，按照0度经线的基准，东经取经度的正值(Longitude)，西经取经度负值(-Longitude)，北纬取90-纬度值(90-Latitude)，南纬取90+纬度值(90+Latitude)，则经过上述处理过后的两点被计为(MLonA, MLatA)和(MLonB, MLatB)。那么根据三角推导，可以得到计算两点距离的如下公式：

C = sin(MLatA)*sin(MLatB)*cos(MLonA-MLonB) + cos(MLatA)*cos(MLatB)

Distance = R*Arccos(C)*Pi/180

这里，R和Distance单位是相同，如果是采用6371.004千米作为半径，那么Distance就是千米为单位，如果要使用其他单位，比如mile，还需要做单位换算，1千米=0.621371192mile
如果仅对经度作正负的处理，而不对纬度作90-Latitude(假设都是北半球，南半球只有澳洲具有应用意义)的处理，那么公式将是：


C = sin(LatA)*sin(LatB) + cos(LatA)*cos(LatB)*cos(MLonA-MLonB)

Distance = R*Arccos(C)*Pi/180

以上通过简单的三角变换就可以推出。
如果三角函数的输入和输出都采用弧度值，那么公式还可以写作：

C = sin(LatA*Pi/180)*sin(LatB*Pi/180) + cos(LatA*Pi/180)*cos(LatB*Pi/180)*cos((MLonA-MLonB)*Pi/180)

Distance = R*Arccos(C)*Pi/180

也就是：

C = sin(LatA/57.2958)*sin(LatB/57.2958) + cos(LatA/57.2958)*cos(LatB/57.2958)*cos((MLonA-MLonB)/57.2958)

Distance = R*Arccos(C) = 6371.004*Arccos(C) kilometer = 0.621371192*6371.004*Arccos(C) mile = 3958.758349716768*Arccos(C) mile

在实际应用当中，一般是通过一个个体的邮政编码来查找该邮政编码对应的地区中心的经纬度，然后再根据这些经纬度来计算彼此的距离，从而估算出某些群体之间的大致距离范围(比如酒店旅客的分布范围-各个旅客的邮政编码对应的经纬度和酒店的经纬度所计算的距离范围-等等)，所以，通过邮政编码查询经纬度这样一个数据库是一个很有用的资源。