動作環境
Xeon E5-2620 v4 (8コア) x 2
32GB RAM
CentOS 6.8 (64bit)
openmpi-1.8.x86_64 とその-devel
mpich.x86_64 3.1-5.el6とその-devel
gcc version 4.4.7 (とgfortran)
NCAR Command Language Version 6.3.0
WRF v3.7.1を使用。
Python 3.6.0 on virtualenv
追記 (2019-01-08)
@knoguchi さんのコメントに記載がありますように、Haversine formula法はGPSウォッチの実測値と少し乖離があるそうです。
より精度の高い計算については 同コメントを参照ください。
情報とv0.1
関連 要検討 > geometry > 緯度経度から距離を求める式の理解 > 極座標から直交座標変換 / 天頂角 / 方位角
緯度経度より距離を(地球の丸さも考えて)求める。のコードをHaversine formula版に変更してみました。
参考 Haversine formulaを使ってローケーションベースでソートする
参考 https://en.wikipedia.org/wiki/Haversine_formula
calc_distance_170331a.py
from math import sin, cos, radians, sqrt, asin
EARTH_RADIUS_km = 6378.137
def dist_on_sphere(pos0, pos1, radius=EARTH_RADIUS_km):
'''
distance based on Haversine formula
Ref: https://en.wikipedia.org/wiki/Haversine_formula
'''
latang1, lngang1 = pos0
latang2, lngang2 = pos1
phi1, phi2 = radians(latang1), radians(latang2)
lam1, lam2 = radians(lngang1), radians(lngang2)
term1 = sin((phi2 - phi1) / 2.0) ** 2
term2 = sin((lam2 - lam1) / 2.0) ** 2
term2 = cos(phi1) * cos(phi2) * term2
wrk = sqrt(term1 + term2)
wrk = 2.0 * radius * asin(wrk)
return wrk
Osaka = 34.702113, 135.494807
Tokyo = 35.681541, 139.767103
London = 51.476853, 0.0
print(dist_on_sphere.__doc__)
print("%.2f km" % dist_on_sphere(Osaka, Tokyo)) # 403.63km
print("%.2f km" % dist_on_sphere(London, Tokyo)) # 9571.22km
実行
$ python calc_distance_170331a.py
distance based on Haversine formula
Ref: https://en.wikipedia.org/wiki/Haversine_formula
403.63 km
9571.22 km
v0.2
radiusの定義を関数内に移した。
キーワード引数で指定すればradiusを変更できるように。
関連 Python > defualt argumentにおいてmutableなobject (例: list)を初期化しない / 定義時に初期化されるため
calc_distance_170331a.py
from math import sin, cos, radians, sqrt, asin
def dist_on_sphere(pos0, pos1, radius=None):
'''
distance based on Haversine formula
Ref: https://en.wikipedia.org/wiki/Haversine_formula
'''
if radius is None:
radius = 6378.137 # km (Earth's radius)
latang1, lngang1 = pos0
latang2, lngang2 = pos1
phi1, phi2 = radians(latang1), radians(latang2)
lam1, lam2 = radians(lngang1), radians(lngang2)
term1 = sin((phi2 - phi1) / 2.0) ** 2
term2 = sin((lam2 - lam1) / 2.0) ** 2
term2 = cos(phi1) * cos(phi2) * term2
wrk = sqrt(term1 + term2)
wrk = 2.0 * radius * asin(wrk)
return wrk
Osaka = 34.702113, 135.494807
Tokyo = 35.681541, 139.767103
London = 51.476853, 0.0
print(dist_on_sphere.__doc__)
print("%.2f km" % dist_on_sphere(Osaka, Tokyo)) # 403.63km
print("%.2f km" % dist_on_sphere(London, Tokyo)) # 9571.22km