计算经纬度两点之间的距离（c#）

(第一部分)

/// <summary>

/// Calculates the distance between to lat/long points and returns the approximate distance in kilometers

/// </summary>

/// <param name="from">Point in long/lat decimal degrees</param>

/// <param name="to">Point in long/lat decimal degrees</param>

/// <returns>Distance in kilometers</returns>

private double CalcDistance(Point from, Point to)

{

    double rad = 6371; //Earth radius in Km

    //Convert to radians

    double p1X = from.X / 180 * Math.PI;

    double p1Y = from.Y / 180 * Math.PI;

    double p2X = to.X / 180 * Math.PI;

    double p2Y = to.Y / 180 * Math.PI;

    

    return Math.Acos(Math.Sin(p1Y) * Math.Sin(p2Y) +

        Math.Cos(p1Y) * Math.Cos(p2Y) * Math.Cos(p2X - p1X)) * rad;

}

 

 

 

 

(第二部分)

 

计算地球上两经纬度点A B间距离

 

在GIS应用中，计算两点之间距离的公式非常重要，这里仅列出几种计算方法。

假设地球是一个标准球体，半径为R,并且假设东经为正，西经为负，北纬为正，南纬为负，

则 A(x,y)的坐标可表示为（R*cosy*cosx,R*cosy*sinx,R*siny） B(a,b)可表示为(R*cosb*cosa,R*cosb*sina,R*sinb)

于是，AB对于球心所张的角的余弦大小为
cosb*cosy*(cosa*cosx+sina*sinx)+sinb*siny

=cosb*cosy*cos(a-x)+sinb*siny
因此AB两点的球面距离为
R*{arccos[cosb*cosy*cos(a-x)+sinb*siny]}

注意几点：

1.      x,y,a,b都是角度，最后结果中给出的arccos因为弧度形式；

2.      所谓的“东经为正，西经为负，北纬为正，南纬为负”是为了计算的方便。比如某点为西经145°，南纬36°，那么计算时可用(-145°,-36°)；

3.      AB对球心所张角的球法实际上是求<OA>和<OB>两向量的夹角K。用公式<OA>*<OB>=|OA|*|OB|*cosK可以得到；

4.      还有对相同点进行处理等。

 

参考资料1给出了计算通过两个点的经纬度计算距离；

原理为：

地球赤道上环绕地球一周走一圈共40075.04公里,而@一圈分成360°,而每1°(度)有60,每一度一秒在赤道上的长度计算如下：

　　40075.04km/360°=111.31955km

　　111.31955km/60=1.8553258km=1855.3m

　　而每一分又有60秒,每一秒就代表1855.3m/60=30.92m

　　任意两点距离计算公式为

　　d＝111.12cos{1/[sinΦAsinΦB十cosΦAcosΦBcos(λB—λA)]}

　　其中A点经度，纬度分别为λA和ΦA，B点的经度、纬度分别为λB和ΦB，d为距离。

c#代码

private const double EARTH_RADIUS = 6378.137; //地球半径

private static double rad(double d)

{

   return d * Math.PI / 180.0;

}

 

public static double GetDistance(double lat1, double lng1, double lat2, double lng2)

{

   double radLat1 = rad(lat1);

   double radLat2 = rad(lat2);

   double a = radLat1 - radLat2;

   double b = rad(lng1) - rad(lng2);

   double 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 * EARTH_RADIUS;

   s = Math.Round(s * 10000) / 10000;

   return s;

}

 

参考资料2给出计算经纬度距离的matlab版本（代码太长，读者可自己链接，这是参考了http://www.ga.gov.au/geodesy/calcs/ 的方法）；

 

    参考资料3给出了从Google Map得到启示的C#版本；

下面就是用C#根据经纬度求两点间距离的函数代码

public static double DistanceOfTwoPoints(double lng1,double lat1,  double lng2, double lat2, GaussSphere gs)
         {            
            double radLat1 = Rad(lat1);
            double radLat2 = Rad(lat2);
            double a = radLat1 - radLat2;
            double b = Rad(lng1) - Rad(lng2);
           double 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 * (gs == GaussSphere.WGS84 ? 6378137.0 : (

gs == GaussSphere.Xian80 ? 6378140.0 : 

6378245.0));
            s = Math.Round(s * 10000) / 10000;
            return s;
        }
        
        private static double Rad(double d)
         {
            return d * Math.PI / 180.0;
        }

    GaussSphere 为自定义枚举类型
    /// <summary>
    /// 高斯投影中所选用的参考椭球
    /// </summary>
    public enum GaussSphere
     {
        Beijing54,
        Xian80,
        WGS84,
    }
 

     参考资料5给出了计算两点经纬度距离的众多方法，给出了计算公式（包括源码）和改进的方法。所有这些公司都是基于地球是球体的假设，这个假设对众多的目的应用已经足够了（实际上地球是一个类似椭球体，用一个球体计算模型最大的误差在0.3%，详见该网页中的笔记部分）。

 

参考资料

1.      某Blog http://blog.csdn.net/yichangxin/archive/2009/02/16/3897553.aspx

2.      经纬度计算距离的matlab版本 http://crust.cn/?p=197

3.      用C#根据经纬度求两点间距离的函数代码 http://www.cnblogs.com/xionglee/articles/1493276.html

4.      权威计算方法 http://www.ga.gov.au/geodesy/calcs/

5.      计算脚本网页 http://www.movable-type.co.uk/scripts/latlong.html