地球上の2点における距離と角度の計算方法

概要

とあるプログラムを書くときに必要になったので。

設定

地球の半径は$R$

地球上の2点$P_i(i=1,2)$は緯度が$\theta_i$で経度が$\varphi_i$である。

このとき、

$$-\pi/2\lt \theta_i\lt \pi/2$$

$$-\pi\lt \varphi_i\lt \pi$$

とする。

原点$O$を地球の中心

z軸を地球の自転軸

つまり、$(0,0,R)$を北極点、$(0,0,-R)$を南極点として、

本初子午線と赤道の交点が$(R,0,0)$

赤道上・東経90°の地点が$(0,R,0)$

となるような座標系を考えてみる。

すると、

$$\overrightarrow{OP_i}=(x_i,y_i,z_i)=(R\cos{\theta_i}\cos{\varphi_i},R\cos{\theta_i}\sin{\varphi_i},R\sin{\theta_i})$$

距離

距離が近いなら、単にベクトルの差のノルムを求めたらそんなに誤差は出ないけど、2点の距離が離れるほど誤差がでかくなるので、きちんと実装しましょう。

地球のような球面だと、最短距離は大円を求めればOK

つまり$O,P_1,P_2$を通る平面で地球を切断したときの扇形$OP_1P_2$の孤の長さを求めればよい。

それは$\angle P_1OP_2$を求めればよいが、それは内積を計算すればいいことがわかる。

つまり、

$$R\cdot \cos^{-1}{\left(\dfrac{\overrightarrow{OP_1}\cdot \overrightarrow{OP_2}}{R^2}\right)}$$

が答え

この式は、P_1,P_2がどこにあっても正しく動くはず

角度

P_1から見たときのP_2の角度を求める。

「P_2から見たときのP_1の角度」ではないことに注意。後述

点$P_1=(x_1,y_1,z_1),x_1^2+y_1^2+z_1^2=R^2$に観測者がいるとき、

真北方向のベクトルは

$$\overrightarrow{N_1}=\left(-\dfrac{x_1z_1}{\sqrt{x_1^2+y_1^2}},-\dfrac{y_1z_1}{\sqrt{x_1^2+y_1^2}},\sqrt{x_1^2+y_1^2}\right)$$

真東方向のベクトルは

$$\overrightarrow{E_1}=\left(\dfrac{y_1R}{\sqrt{x_1^2+y_1^2}},-\dfrac{x_1R}{\sqrt{x_1^2+y_1^2}},0\right)$$

$\overrightarrow{E_1}$ベクトルと$\overrightarrow{N_1}$ベクトルは共にノルムが$R$であることに注意する。

（注：$P_1$が北極点または南極点である時、$\sqrt{x_1^2+y_1^2}=0$であるため分母に0が入って正常に計算できない。北極点から見たときはあらゆる方角が南、南極点から見たときはあらゆる方角が北であることに注意せよ）

このとき

$$Y=\overrightarrow{N_1}\cdot \overrightarrow{P_1P_2}$$

$$X=\overrightarrow{E_1}\cdot \overrightarrow{P_1P_2}$$

とすると、$(0,0)$から$(X,Y)$への方角と等しくなる

よって、ATAN2関数に$(X,Y)$を代入すれば答えが出る。

ここで、$P_2$がちょうど$P_1$の対蹠点にある場合については、$\overrightarrow{P_1P_2}$が$\overrightarrow{N_1},\overrightarrow{E_1}$と直交するので、$(X,Y)=(0,0)$となる。この調子でプログラムを組むと、ATAN2関数の仕様に沿った答えが出力される。ところで、この場合にはどの方向に進んでも$P_2$に到達できるので、当然ATAN2(0,0)の方向に進んでも目的地に到達できるわけである。よって、答えとしては間違ってはいないことになる。

結局、このプログラムを組んだ場合、$P_1$が北極点or南極点でない限り正しい答えが出る。