Python 経度・緯度で与えられた２点間距離計算
はじめに
この記事は、Qiita に投稿した以下の記事の再アレンジ版です。
経度・緯度で与えられた２点間距離を計算するプログラムです。 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点 と B点
の測地線長 (
) を求める。
計算手順
（１）測地緯度を化成緯度に変換
（２）球面上の距離 (Spherical Distance:
)の計算
（３）距離補正量 (Lambert-Andoyer Correction)
（４）測地線長(Geodetic Distance)
プログラム例
- 経度・緯度の情報は，小数点方式で与えます．度・分・秒形式のデータの場合には，度＋分➗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
