已知经纬度坐标求亮点距离

此文转载

例题:

已知巴黎的地理位置为 东经2°20'14'',北纬 48°50'11'，华盛顿的地理位置为 西经 22°03'56'',北纬 38°55'17''
求两城市之间的距离。

解答:
用低精度和高精度两种方法解答，两种方法的区别在于，低精度将地球看作球形，高精度将地球看作椭球。
首先将经纬度转换成十进制的度。
定义函数deg2dec，用来完成这一转换

deg2dec <-

function(h,m,s){

    if(h < 0)

    {m = - m

     s = -s

    }

    res = h + m/60 + s/3600

    return(res)

}

 

##巴黎

L1 = deg2dec(-2,20,14)

phi1 = deg2dec(48, 50, 11)

##华盛顿

L2 = deg2dec(77,03,56)

phi2 = deg2dec(38,55,17)

 

 

## 第一种低精度方案，

## 由于地球可以的扁率较小，可以近似看作一个球体，利用球面三角学公式可以快速的求出两点之间的距离。写成计算机程序，如下所示:

 

#低精度公式

geodist <-

function(L1, phi1, L2, phi2)

{

   Ri = 6371

   d = ( Ri * acos(sin(phi1*(pi/180))*sin(phi2*(pi/180)) + cos(phi1*(pi/180))*cos(phi2*(pi/180))*cos((L1 - L2)*(pi/180))))

   res <- round(d, 3)

   return(res)

}

 

#低精度

geodist(L1, phi1, L2, phi2)

答案：

华盛顿和巴黎之间的距离是 6165.597km

 

##第二种 高精度方案，高精度的方案的误差低于+_50m.

## 采用较为精密的公式，考虑到地球是一个椭球,a是地球的长轴，f是扁率，用如下公式计算

 

 

hageodist <-

function(L1, phi1, L2, phi2){

    a = 6378.14

    f = 1/298.257

    F = (phi1+phi2)/2

    G = (phi1 - phi2)/2

    ramda <- (L1 - L2)/2

   

    S = (sin(G*pi/180)^2)*(cos(ramda*pi/180)^2) + (cos(F*pi/180)^2)*(sin(ramda*pi/180)^2)

    C= (cos(G*pi/180)^2)*(cos(ramda*pi/180)^2) + (sin(F*pi/180)^2)*(sin(ramda*pi/180)^2)

   

    omega = atan(sqrt(S/C))

    R = sqrt(S*C)/omega

    D = 2*omega*a

   

    H1 = (3*R-1)/(2*C)

    H2 = (3*R+1)/(2*S)

   

    res = D*(1 + f*H1*(sin(F*pi/180)^2)*(cos(G*pi/180)^2) - f*H2*(cos(F*pi/180)^2)*(sin(G*pi/180)^2))

 return(round(res,3))

}

 

hageodist(L1, phi1, L2, phi2)

答案：

华盛顿和巴黎之间的距离是6181.628km.

在计算这一距离时，高精度和低精度算法给出的结果之间差距达16km。可见在天文计算中，低精度是远远不能满足要求的。

 

 

例题和公式参考 Jean Meeus 的 Astronomical Algorithms