Python
距離計算

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

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

緯度・経度で与えられた2点間距離を計算するプログラム。
Lambert-Andoyerの方法を用いています。
参考サイトは以下の通り。

測地線長の計算(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} \qquad \text{$F$:偏平率 (Flattening of the Earth)}
\end{equation*}
$$

(4)測地線長(Geodetic Distance)

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

プログラム例

  • 経度・緯度の情報は,小数点方式で与えます.度・分・秒形式のデータの場合には,度+分➗60+秒➗3600 です.
  • 角度は,全てラジアンに換算します.
from math import *


def cal_phi(ra,rb,lat):
    return atan(rb/ra*tan(lat))


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


def main():
    for iii in (0,1):
        if iii==0:
            lat_a= 3.117; lon_a=101.550; loc_a='Kuala Lumpur'
            lat_b=35.690; lon_b=139.760; loc_B='Tokyo'
        if iii==1:
            lat_a= 3.117; lon_a=101.550; loc_a='Kuala Lumpur'
            lat_b=-6.183; lon_b=106.833; loc_B='Jakarta'

        rho=cal_rho(lat_a,lon_a,lat_b,lon_b)
        print('(lat_a, lon_a)=({0:8.3f},{1:8.3f}): {2:16s}'.format(lat_a,lon_a,loc_a))
        print('(lat_b, lon_b)=({0:8.3f},{1:8.3f}): {2:16s}'.format(lat_b,lon_b,loc_B))
        print('Distance={0:10.3f} km'.format(rho))
        print()


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

実行結果

(lat_a, lon_a)=(   3.117, 101.550): Kuala Lumpur    
(lat_b, lon_b)=(  35.690, 139.760): Tokyo           
Distance=  5332.840 km

(lat_a, lon_a)=(   3.117, 101.550): Kuala Lumpur    
(lat_b, lon_b)=(  -6.183, 106.833): Jakarta         
Distance=  1184.232 km