はじめに

地球上の特定の2点の緯度経度座標が与えられ、この2点間の距離を求めたいとします。
ただし考慮するのは緯度経度のみとして、標高は考慮しません。

一般的に地球は完全な球ではなく赤道方向に扁平な回転楕円体として知られています。
よってその表面の2点間の距離を求めるのはなかなか複雑な計算が求められます。

今回はヒュベニの公式と呼ばれるものを採用します。

ヒュベニの公式

かなりの長距離や高緯度地域での使用では誤差が大きく出てくるが、日本国土においての利用程度なら十分な精度で距離が得られそうなので今回採用。

出典

どうにも出典がはっきりしませんが、「カシミール3D」という登山者向けPCソフトの距離計算式に用いられたことで日本国内での知名度を得たようである。[1]
その際の出典は [2] とされるようだ。

また、次のような記述も見られる。([3])

通常，日本国内で「ヒュベニの公式」や「ヒュベニの距離計算式」などと呼ばれているものは，ヒュベニ（Karl Hubeny）が考案した「ガウスの平均緯度式の改良式」の第1項のみを使用したものをいいます。

さらに、wikipedia が出典になってしまうが、英語版 wikipedia の Geographical distance [4] の中に次のような記述がみられるので、本邦においてヒュベニの公式と呼ばれるものはこれのことを指しているのだろう。
これによると、ヒュベニの公式の出典は [5] ということになる。

Karl Hubeny got the expanded series of Gauss mid-latitude one represented as the correction to flat-surface one.

また、この第1項のみを使用しているという点から、一般にみられるヒュベニの公式を簡易式や簡略式と位置し、ガウスの平均緯度式の拡張を正式なヒュベニの公式としている例もみられるが(([3]),[6])、短距離なら十分な精度が出せるようなので、今回はこの場合のいわゆる簡易式で実装することにする。

式

緯度経度を用いた3つの距離計算方法 [1] に示される式を引用する。

\begin{align}

s & = \sqrt{(M \Delta \phi)^2 + (N \cos \phi \Delta \lambda)^2} \\

\phi & : 2点の平均緯度 \\

\Delta \phi & : 2点の緯度差(ラジアン) \\

\Delta \lambda &: 2点の経度差(ラジアン) \\

M & : 子午線曲率半径 \\
& = a(1-e^2) / \sqrt{(1-e^2 \sin^2 \phi)^3} \\

N & : 卯酉線曲率半径 \\
& = a / \sqrt{1-e^2 \sin^2 \phi} \\

a & : 回転楕円体の長半径(赤道半径) \\

b & : 回転楕円体の短半径(極半径) \\

e & : 離心率 = \sqrt{1-(b/a)^2}

\end{align}

この時、赤道半径・離心率は GRS80 [7] からそれぞれ以下の固定値としておく。
極半径は離心率の算出以外に使用しないので省略。

\begin{align}

a & : 回転楕円体の長半径(赤道半径) = 6378137 (m) \\
e & : 離心率 = 0.081819191042815791

\end{align}

実装

ある2点の緯度と経度が得られたときに、その距離を計算するコードを Ruby の method で実装する

def calc_hubeny_distance(lat_x, lng_x, lat_y, lng_y)
  # ラジアンに換算
  lat_x_rad = lat_x * Math::PI / 180
  lng_x_rad = lng_x * Math::PI / 180
  lat_y_rad = lat_y * Math::PI / 180
  lng_y_rad = lng_y * Math::PI / 180

  lat_d = lat_x_rad - lat_y_rad # 緯度差
  lng_d = lng_x_rad - lng_y_rad # 経度差

  lat_avg = (lat_x_rad + lat_y_rad) / 2.0 # 緯度平均

  e = 0.081819191042815791 # e: 地球の離心率 (GRS80に従う)

  equatorial_radius = 6378137.0 # a: 赤道半径(m) (GRS80に従う)

  w = Math.sqrt(1.0 - (e ** 2) * Math.sin(lat_avg) ** 2) # 子午線・卯酉線曲率半径の分母

  m = equatorial_radius * (1.0 - e ** 2) / w ** 3 # 子午線曲率半径
  n = equatorial_radius / w # 卯酉線曲率曲率半径

  t1 = lat_d * m
  t2 = lng_d * n * Math.cos(lat_avg)

  Math.sqrt(t1 ** 2 + t2 ** 2)
end

横着をして method にしたので固定値とした離心率と赤道半径が変数っぽいが気にしないことにする。

さいごに

最後に、東京駅 ~ 大阪駅間の距離を算出して結果のほどを確かめてみる。
今回は国土地理院の距離と方位角の計算のページの結果を正とする。

# 東京駅
point_x = {lat: 35.6809591, lng: 139.7673068}
# 大阪駅
point_y = {lat: 34.7022887, lng: 135.4953509}

delta = calc_hubeny_distance(
  point_x[:lat],
  point_x[:lng],
  point_y[:lat],
  point_y[:lng]
)

p "#{delta} m"
# => "403934.14283672185 m"

また、先述の国土地理院のページで同じ座標間の距離を求めたところ 403,893.941(m) となった。

日本くらいの中程度の緯度地域では 400 km 程度の距離で 40 m 程の誤差が出る可能性があるようだが、簡単な実装にしては十分な精度が出ているだろう。

参考・引用

[1] 三浦　英俊. 緯度経度を用いた3つの距離計算方法. オペレーションズ・リサーチ 60(12) 701-705, 2015
[2] 高崎正義，『地図学』，朝倉書店，1988.
[3] 歩測計 Pace Ruler - ASTI アマノ技研. FAQ（よくある質問と回答）
[4] wikipedia - Geographical distance
[5] Hubeny, K. (1954). Entwicklung der Gauss'schen Mittelbreitenformeln, Österreichische Zeitschrift für Vermessungswesen.
[6] mapli.net - Distance 地図上の2点間の最短距離と方位角を楕円面上で精密に求める
[7] Wikipedia - GRS80