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

(2018.09.10) プログラムを．if __name__ == '__main__': main()方式で書き直しました．
(2019.05.19) numpy対応にしました、また文書の記載ミスを修正しました。

はじめに

経度・緯度で与えられた２点間距離を計算するプログラムです。
Lambert-Andoyerの方法を用いています。
データ入力の並びは、GMTのプロットに合わせて、「経度・緯度」の順にしています。

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

• Lambert-Andoyerの方法：http://www2.nc-toyama.ac.jp/~mkawai/lecture/sailing/geodetic/geosail.html
• おもしろい実験：http://petit-noise.net/blog/%E7%B7%AF%E5%BA%A6%E7%B5%8C%E5%BA%A6%E3%81%8B%E3%82%892%E7%82%B9%E9%96%93%E3%81%AE%E8%B7%9D%E9%9B%A2%E3%82%92%E6%B1%82%E3%82%81%E3%82%8B/
• 赤道半径と極半径：https://www.jamstec.go.jp/j/kids/qa/02_1.html

測地線長の計算（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}

計算手順

（１）測地緯度を化成緯度に変換

\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}

（２）球面上の距離 (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}

（３）距離補正量 (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}

（４）測地線長(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

以　上