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

最新推荐文章于 2024-02-03 22:41:13 发布

jasonCSH

最新推荐文章于 2024-02-03 22:41:13 发布

阅读量6.6k

点赞数

文章标签： matlab c# google 脚本 blog 360

在 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,
    }