[Ruby]緯度と経度から2点間の距離を求める(ヒュベニの公式)

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

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

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

# ヒュベニの公式

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

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

短距離なら十分な精度が出せるようなので、今回はこの式で実装することにする。

## 式

緯度経度を用いた3つの距離計算方法 [[1]](https://orsj.org/wp-content/corsj/or60-12/or60_12_701.pdf) に示される式を引用する。

```math
\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](https://ja.wikipedia.org/wiki/GRS80)] からそれぞれ以下の固定値としておく。
極半径は離心率の算出以外に使用しないので省略。

```math
\begin{align}

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

\end{align}
```

# 実装

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

```ruby
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 にしたので固定値とした離心率と赤道半径が変数っぽいが気にしないことにする。

# さいごに
最後に、東京駅 ~ 大阪駅間の距離を算出して結果のほどを確かめてみる。
今回は国土地理院の [距離と方位角の計算](https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2stf.html) のページの結果を正とする。

```ruby
# 東京駅
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](https://orsj.org/wp-content/corsj/or60-12/or60_12_701.pdf)
[2] 高崎正義，『地図学』，朝倉書店，1988.
~~[3] [歩測計 Pace Ruler - ASTI アマノ技研. FAQ（よくある質問と回答）](https://amano-tec.com/apps/paceruler.html#jp_faq)~~
~~[4] [wikipedia - Geographical distance](https://en.wikipedia.org/wiki/Geographical_distance)~~
~~[5] Hubeny, K. (1954). Entwicklung der Gauss'schen Mittelbreitenformeln, Österreichische Zeitschrift für Vermessungswesen.~~
~~[6] [mapli.net - Distance 地図上の2点間の最短距離と方位角を楕円面上で精密に求める](https://www.mapli.net/distance/)~~
[7] [Wikipedia - GRS80](https://ja.wikipedia.org/wiki/GRS80)