概要
地球上の2地点間の距離を求めることってたまにありますよね．ただ，その2地点の情報は経度緯度で与えられることが多く，単純には距離を求めることはできません．
経度緯度から距離を求める代表的な方法として，測地線航海算法（Geodesic Sailing）というものがあります．今回はPythonを使って，測地線航海算法による2地点間の距離を計算してみましょう．
測地線航海算法とは
こちらの富山高等専門学校の解説サイトを参考に解説します．
公式
結論から言うと，（緯度，経度）がそれぞれ$(l_1,L_1),(l_2,L_2)$である2地点間の距離$\rho$は，次の式で求められます．ただし$A$は地球赤道半径，$B$は地球極半径です．
\begin{align}
\phi_1 &= \tan^{-1}\left(\frac{B}{A}\tan l_1\right),\phi_2 = \tan^{-1}\left(\frac{B}{A}\tan l_2\right)\tag{1}\\
X &= \cos^{-1}\{\sin\phi_1 \sin\phi_2 + \cos\phi_1 \cos\phi_2 \cos(L_1-L_2)\}\tag{2}\\
F &= \frac{A-B}{A}\tag{3}\\
\Delta\rho &= \frac{F}{8} \left\{ (\sin X - X)\left(\frac{\sin\phi_1+\sin\phi_2}{\cos\frac{X}{2}}\right)^2 - (\sin X + X)\left(\frac{\sin\phi_1-\sin\phi_2}{\sin\frac{X}{2}}\right)^2 \right\}\tag{4}\\
\rho &= A(X+\Delta\rho)
\end{align}
解説（概略）
求める距離$\rho$は測地線長（Geodesic Distance）と呼ばれ，回転楕円体（≒地球）上の2地点間の最短距離を指します．
1. 化成緯度への変換
まず式$(1)$では，我々がふつう経度と呼ぶ測地緯度（地理緯度）から化成緯度というものに変換しています．
ある地点$X$の測地緯度とは，下図のように$X$における楕円体の法線を引いたときの角度$l_A$のことです．
画像引用：http://motochan.sakura.ne.jp/public_html/GISContents/16.htm
一方，ある地点$X$の化成緯度とは，下図のように楕円体上$X$を球面上の$X'$に移し，$X'$における法線を引いた時の角度$l_D$のことを指します．
この$l_A$と$l_D$の関係が式$(1)$で表せるというわけです．
画像引用：http://motochan.sakura.ne.jp/public_html/GISContents/16.htm
2. 球面上の距離の計算
次に式$(2)$では，球面に移した点の間の距離$X$（ラジアン）を計算しています．
簡単のため$L_1=L_2$の場合を考えると，球面に移した2点の間の距離$X$は単に
$$X=\phi_1−\phi_2$$
となります（ここで言う距離の単位は角度（rad or deg））．
一方，式$(2)$で$L1=L2$とすると，
\begin{align}
X &= \cos^{-1}\{\sin\phi_1 \sin\phi_2 + \cos\phi_1 \cos\phi_2\}\\
&= \cos^{-1}\cos(\phi_1 - \phi_2)\\
&= \phi_1 - \phi_2
\end{align}
となり，合致しています．
$L_1\ne L_2$の場合はそう単純にはいかないため，$\cos(L_1-L_2)$がついているようです．
3. 扁平率
式$(3)$の$F$は扁平率と呼ばれ，長半径と短半径の差を長半径で割った値で定義されます．
なので，扁平率は大きければ大きいほど潰れた形をしています．
地球の場合，扁平率は約$\frac{1}{300}$なのでほぼ球に近いと言えます．
4. Lambert-Andoyer補正
式$(4)$では，地球が球面ではなく楕円体であることを補正するための計算を行っています．
もし地球が球面なら，扁平率$F=0$となるので，式$(4)$は$\Delta\rho=0$となります．
5. 測地線距離
あとは，2.で求めた$X$に補正項$\Delta\rho$を足した$X+\Delta\rho$に半径$A$をかけてやれば1，求める距離$\rho$が得られます．
Pythonで書いてみる
Pythonで測地線航海算法により距離を求めてみましょう．
import numpy as np
# 世界測地系(GRS80)
A = 6378137.0
F = 1/298.257222101
B = A*(1-F)
# 日本測地系(Bessel)
# A = 6377397.155
# F = 1/299.152813
# B = A*(1-F)
# Google Maps
# A = 6371008
# B = A
# F = (A-B)/A
def dist(p1, p2):
global A, B, F
# 緯度経度をラジアンに変換
lat1 = np.radians(p1[0])
lon1 = np.radians(p1[1])
lat2 = np.radians(p2[0])
lon2 = np.radians(p2[1])
# 化成緯度に変換
phi1 = np.arctan2(B * np.tan(lat1), A)
phi2 = np.arctan2(B * np.tan(lat2), A)
# 球面上の距離
X = np.arccos(np.sin(phi1) * np.sin(phi2) + np.cos(phi1) * np.cos(phi2) * np.cos(lon2 - lon1))
# Lambert-Andoyer補正
drho = F/8 * ((np.sin(X) - X) * (np.sin(phi1) + np.sin(phi2))**2 / np.cos(X/2)**2 - (np.sin(X) + X) * (np.sin(phi1) - np.sin(phi2))**2 / np.sin(X/2)**2)
# 距離
rho = A * (X + drho)
return rho
試しに，北海道庁から沖縄県庁までの距離を計算してみましょう．
# 都道府県庁の緯度経度は https://www.gsi.go.jp/common/000230936.pdf を参考
if __name__ == '__main__':
p1 = (43+3/60+52/3600, 141+20/60+49/3600) # 北海道庁
p2 = (26+12/60+45/3600, 127+40/60+51/3600) # 沖縄県庁
print(dist(p1,p2))
結果は以下の通り．
ほぼ正確に計算できていることがわかります．
| 距離[m] | |
|---|---|
| 自分の実装 | 2243872.655854546 |
| 国土地理院の測量計算 | 2243875.695 |
ソースコード
今回は大したコードは書いていませんが，Qiitaアドベントカレンダーで使ったコードはこちらのレポジトリにまとめています．
本記事のコードは04-distanceの中に入っています．
-
円弧の長さの公式：$l=r\theta$を思い出してください． ↩