geometry > 2017-03-31: 2つの緯度経度から距離を求める > Haversine formula使用 | 2019-01-08 Vincenty法

動作環境

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