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

緯度・経度で与えられた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

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()

実行結果

(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
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account log in.