計算地球上兩經緯度點 A B 間距離

現在的位置: 首頁 > 綜合 > 正文

RSS

計算地球上兩經緯度點 A B 間距離

2018年04月18日 ⁄ 綜合 ⁄ 共 3103字 ⁄ 字型大小小中大 ⁄ 評論關閉

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