Edited at

Python 経度・緯度で与えられた2点間距離計算

(2018.09.10) プログラムを.if __name__ == '__main__': main()方式で書き直しました.

(2019.05.19) numpy対応にしました、また文書の記載ミスを修正しました。


はじめに

経度・緯度で与えられた2点間距離を計算するプログラムです。

Lambert-Andoyerの方法を用いています。

データ入力の並びは、GMTのプロットに合わせて、「経度・緯度」の順にしています。

参考サイトは以下の通り。


測地線長の計算(Lambert-Andoyer Method for Computation of Geodetic Distance)


設問

A点 $(lon_A,lat_A)$ と B点 $(lon_B,lat_B)$ の測地線長 ($\rho$) を求める。

\begin{align}

&\text{ここに}& &(lon_A,lat_A)& &\text{:A点の測地経度・測地緯度} \\
& & &(lon_B,lat_B)& &\text{:B点の測地経度・測地緯度}
\end{align}


計算手順


(1)測地緯度を化成緯度に変換

\begin{equation}

\phi=\tan^{-1}\left(\cfrac{r_b}{r_a} \cdot \tan (lat) \right)
\end{equation}

\begin{align}

&\phi & &\text{:化成緯度 (Parametric Latitude)}\\
&lat & &\text{:測地緯度 (Geodetic Latitude)}\\
&r_a & &\text{:地球赤道半径 (Semi-major Axis of the Earth)}\\
&r_b & &\text{:地球極半径 (Semi-minor Axis of the Earth)}
\end{align}


(2)球面上の距離 (Spherical Distance:X)の計算

\begin{equation}

X=\cos^{-1}\left(\sin\phi_A \cdot \sin\phi_B+\cos\phi_A \cdot \cos\phi_B \cdot \cos(lon_A-lon_B)\right)
\end{equation}

\begin{align}

&\phi_A & &\text{:A点の化成緯度}\\
&\phi_B & &\text{:B点の化成緯度}
\end{align}


(3)距離補正量 (Lambert-Andoyer Correction)

\begin{equation}

\Delta\rho=\cfrac{F}{8}\cdot\left\{\cfrac{[(\sin X-X)(\sin\phi_A+\sin\phi_B)]^2}{\cos^2 (X/2)}
-\cfrac{[(\sin X+X)(\sin\phi_A-\sin\phi_B)]^2}{\sin^2 (X/2)}\right\}
\end{equation}

\begin{equation}

F=\cfrac{r_a-r_b}{r_a} \quad \text{:偏平率 (Flattening of the Earth)}
\end{equation}


(4)測地線長(Geodetic Distance)

\begin{equation}

\rho=r_a\cdot(X+\Delta\rho) \quad \text{:測地線長}
\end{equation}


プログラム例


  • 経度・緯度の情報は,小数点方式で与えます.度・分・秒形式のデータの場合には,度+分➗60+秒➗3600 です.

  • 角度は,全てラジアンに換算します.

  • プログラム例では、クアラルンプールから各都市までの距離を計算しています。

  • データはリスト形式で [経度, 緯度, 都市名] という順に与えています。

import numpy as np

def cal_rho(lon_a,lat_a,lon_b,lat_b):
ra=6378.140 # equatorial radius (km)
rb=6356.755 # polar radius (km)
F=(ra-rb)/ra # flattening of the earth
rad_lat_a=np.radians(lat_a)
rad_lon_a=np.radians(lon_a)
rad_lat_b=np.radians(lat_b)
rad_lon_b=np.radians(lon_b)
pa=np.arctan(rb/ra*np.tan(rad_lat_a))
pb=np.arctan(rb/ra*np.tan(rad_lat_b))
xx=np.arccos(np.sin(pa)*np.sin(pb)+np.cos(pa)*np.cos(pb)*np.cos(rad_lon_a-rad_lon_b))
c1=(np.sin(xx)-xx)*(np.sin(pa)+np.sin(pb))**2/np.cos(xx/2)**2
c2=(np.sin(xx)+xx)*(np.sin(pa)-np.sin(pb))**2/np.sin(xx/2)**2
dr=F/8*(c1-c2)
rho=ra*(xx+dr)
return rho

def main():
no=[[101.550, 3.117, 'Kuala Lumpur'],
[139.760, 35.690, 'Tokyo'],
[106.833, -6.183, 'Jakarta'],
[100.567, 13.733, 'Bangkok'],
[102.567, 17.950, 'Vientiane'],
[105.800, 21.017, 'Ha Noi'],
[126.967, 37.567, 'Seoul'],
[116.283, 39.933, 'Beijing'],
[120.967, 14.583, 'Maynila']]
n=8
loc_a=[no[0][2]]*n
lon_a=np.full(n,no[0][0])
lat_a=np.full(n,no[0][1])
loc_b= [no[1][2],no[2][2],no[3][2],no[4][2],no[5][2],no[6][2],no[7][2],no[8][2]]
lon_b=np.array([no[1][0],no[2][0],no[3][0],no[4][0],no[5][0],no[6][0],no[7][0],no[8][0]])
lat_b=np.array([no[1][1],no[2][1],no[3][1],no[4][1],no[5][1],no[6][1],no[7][1],no[8][1]])

rho=cal_rho(lon_a,lat_a,lon_b,lat_b)
print('{0:<16s} {1:<16s} {2:>8s} {3:>8s} {4:>8s} {5:>8s} {6:>12s}'\
.format('loc_a','loc_b','lon_a','lat_a','lon_b','lat_b','Dist.(km)'))
for i in range(len(loc_a)):
print('{0:16s} {1:16s} {2:8.3f} {3:8.3f} {4:8.3f} {5:8.3f} {6:12.3f}'\
.format(loc_a[i],loc_b[i],lon_a[i],lat_a[i],lon_b[i],lat_b[i],rho[i]))

#==============
# Execution
#==============
if __name__ == '__main__': main()


実行結果

loc_a            loc_b               lon_a    lat_a    lon_b    lat_b    Dist.(km)

Kuala Lumpur Tokyo 101.550 3.117 139.760 35.690 5332.840
Kuala Lumpur Jakarta 101.550 3.117 106.833 -6.183 1184.232
Kuala Lumpur Bangkok 101.550 3.117 100.567 13.733 1179.107
Kuala Lumpur Vientiane 101.550 3.117 102.567 17.950 1644.533
Kuala Lumpur Ha Noi 101.550 3.117 105.800 21.017 2033.151
Kuala Lumpur Seoul 101.550 3.117 126.967 37.567 4612.783
Kuala Lumpur Beijing 101.550 3.117 116.283 39.933 4339.821
Kuala Lumpur Maynila 101.550 3.117 120.967 14.583 2480.551

以 上