今回は懐かしの地理の知識と数学の融合である。飛行機で欧米に行ったことがある人なら見覚えがあると思うが、世界地図上で飛行機の航路が北半球なら上に膨らむように取ることが多い。
国際情勢上やジェット気流の影響は考えないものとして、最短航路を描画するにはどんなアルゴリズムを作ればいいだろうか？興味本位でやってみよう。

書いてあることはwikiの焼き直しに近いが、数学の力をもってすれば日本語の少ない要件を覚えるだけで、数式系は暗記せずともマッピングのルールを導くことができるのである。数学というツールは恐るべし。

なお、wikiに書いてある数式は緯度や経度の基準点をどこにおくかでも変わってしまう点もあるので、自力実装するならば、数式は暗記や参照せずに自力で考え方の理解または導出して使った方が安全だといえる。

球面を平面内に描写する図法

中学校地理の授業のおさらいである。地球儀を一枚の平面上に地図として表現する図法には、正射図法、メルカトル図法、サンソン図法、正距方位図法などがある(他にモルワイデ図法が使われることが多いが今回は
説明を省略する)。

今回やりたいのは、

(世界地図引用するのが面倒なので)任意の画像を歪み最小限で球面上にマッピングする。
マッピングした球面上で今回のお題である、2点間のメルカトル図法上での最短航路を描画する
の順に行う。これに便利なのが先ほどの図法ということだ。

各図法の要件と定義式

球座標

今回は計算の便宜上、北極を0度とする一般的な極座標で表現する。
(こうした場合、wikiに書いてある計算式はそのまま使うことはできなくなる点は注意であるが)

$$x=\sin\theta\cos\phi$$
$$y=\sin\theta\sin\phi$$
$$z=\cos\theta$$

正射図法

無限遠方から球を眺めたときの像を描写する。
北極を0度とする

とくに、いまy軸方向から球を眺めた場合、カメラ座標(u,v)は、カメラパラメータを適当に設計すれば

$$u=x=\sin\theta\cos\phi$$
$$v=z=\cos\theta$$

で表現できる。

サンソン図法

球面のパターンを経線方向に輪切りにして、ある緯度を中心に揃えて縦に並べた図法である。地図上の点(u,v)は、縦軸は緯度のみに依存し、横軸は経度方向を輪切りの周長に沿って描画する。

$$u=\phi\sin\theta$$
$$v=\theta$$

このとき、任意の微小面積(du,dv)は球面上で面積が等しくなるように対応する。このような図法を等積図法とよぶ。

メルカトル図法

以下の2点を要件とする座標である。

緯度と経度をどこでも直角に交差するようにする
かつ出発点の角度は平行移動しても維持するようにする

2の特性をもつ図法を等角図法と呼ぶ。
これに対応する変換式を導こう。
1の特性を満たす座標(u,v)は、任意の関数$f,g$を用いて

$$u=f(\phi)$$
$$v=g(\theta)$$

で表現できる。つまり、$u,v$の座標が$\phi,\theta$いずれか一方ずつ依存していれば直交性は保証される。
あとは等角性を保つように$f,g$を定めればよい。経度のスケールを固定するため$f(\phi)=\phi$
として$g$を求めよう。

もし経度方向を全緯度において正しいスケールで描画するには輪切りの周長を考慮しなければならなかったが、要件1によって直行性を保つためにスケールの不変性を捨てている。その倍率は緯度方向に依存し、
$$u=\phi\sin\theta\cdot s(\theta)=\phi$$
$$\therefore s(\theta)=\frac{1}{\sin\theta}$$
である。
要件2の等角性を保つには、縦軸と横軸の倍率は同じでなければならない。つまりある点$(u,v)$における角度微小変化$d\phi,d\theta$に対する変化量$(du,dv)$は、

$$du=\sin\theta s(\theta)d\phi=d\phi$$
$$dv=s(\theta)d\theta$$

である必要がある。定数を適当に設定すれば$u=\phi$は自明である。$v$は、

$$v=\int s(\theta)d\theta$$
となる。これを解こう。

$$v=\int \frac{1}{\sin\theta}d\theta$$

$$=\int \frac{\sin\theta}{1-\cos^2\theta}d\theta=\int \frac{1}{t^2-1}dt=\frac{1}{2}\int \frac{1}{t-1}-\frac{1}{t+1}dt=\ln\sqrt\frac{1-\cos\theta}{1+\cos\theta}$$
$$\therefore v=\ln\tan\frac{\theta}{2}-v_0\quad (0<\theta<\pi)$$

ここで、$v_0$は$\theta=\frac{\pi}{2}$のときに$v=0$となるように決めよう。

$$v=\ln\tan\frac{\theta}{2}-\ln\tan\frac{\pi}{4}=\ln\tan\frac{\theta}{2}$$

なお、$\theta$の代わりに赤道を0とする緯度を用いる場合は
$$\theta'=\frac{\pi}{2}+\theta$$
と置き換えれば、
$$v=\ln\tan\left(\frac{\pi}{4}+\frac{\theta'}{2}\right)$$
となり、wikiに記載の表式になる。

正距方位図法

これは、球面上の1点を中心点に据えた際に任意の点への方位と距離が等しくなるような図法である。
この要件を満たすには、以下の変換を行えばよい。

中心点が極になるように座標変換(回転)を行う
新しい座標上で同じ緯度の点を同心円上に経度を角度と対応づけて展開する。これで正しい方位関係となる
動径方向のスケールは緯度の定数倍で規定する

北極を中心とする場合

まずは、北極を中心点とした描画を行う際の座標変換をやってみよう。1は不要なので2から立式を始める。

$$u=r\cos\phi$$
$$v=r\sin\phi$$
$$r=k\theta$$

$k$は図の縮尺であり、ここまでの議論の延長で1に固定しよう。つまり
$$u=\theta\cos\phi$$
$$v=\theta\sin\phi$$
が対応することになる。

任意の点を中心におく場合

あとは、1の座標変換で任意の点が極になるような座標に変換する式を求めればよい。
座標変換は三次元座標(x,y,z)ベースで考えた方が簡単である。目標は、$\theta,\phi$の角度に対応する点(x,y,z)が座標回転後に(0,0,1)となるように変換する式である。
このような回転方法は無数にあるので、「着目点と同じ経度の点は変換後に経度0とするように変換」を条件に加えよう。
これを実現する回転は、以下の手順で実現できる。

z軸周り回転を行い、着目点の経度が0になるように合わせる(着目点は変換後にx-z平面内に存在する)
y軸周り回転を行って、着目点をz軸に合わせる

座標回転は過去記事(主成分分析の記事)で議論したとおり、それぞれ回転行列で考えることができる。

$$U_z=\exp(J_z\theta_z)=\begin{pmatrix}\cos\theta_z & \sin\theta_z & 0\
-\sin\theta_z & \cos\theta_z & 0\
0 & 0 & 1
\end{pmatrix}$$

回転行列の指数表記も主成分分析の記事参照である。

$$U_x=\exp(J_x\theta_x),U_y=\exp(J_y\theta_y)$$は$U_z=\exp(J_z\theta_z)$を対角線方向に回すように置き換えることで得られる。

このようにして得られた回転行列$U_y,U_z$を用いて、

$$X'=U_yU_zX=\exp(J_y\theta_y)\exp(J_z\theta_z)X=UX$$

とかける。ただし、$$X=(x\quad y\quad z)^T$$である。

$\theta_y,\theta_z$は着目点(球座標上で$\theta,\phi$)が変換後に$(0,0,1)$となるように定めるので、

$$\theta_z=\phi,\theta_y=\theta$$

とすればよいことがわかる(行列積をとれば検算できる)。

回転後の座標$(x',y',z')$が

$$x'=\sin\theta\cos\phi$$
$$y'=\sin\theta\sin\phi$$
$$z'=\cos\theta$$

を満たすように$\rho,\psi$を決めればよい。

$$\rho=\arccos z'\quad (0\leq\rho\leq\pi)$$
$$\psi=\arccos \frac{x'}{\sin\rho+\epsilon}\mathrm{sign}(y')$$

ただし、$\epsilon$はゼロ除算を防ぐための微小定数である。

求めた角度$\rho,\psi$を用いて

$$u=\rho\cos\psi$$
$$v=\rho\sin\psi$$

でマッピングすれば、正距方位図法が得られる。
正距方位図法上では、中心点から任意の点への最短コースは直線で得られる(大圏コースと呼ぶ)。

最短コースをメルカトル図法上で表す

最短コースはメルカトル図法上では一般に曲線で結ぶことになると想定されるので、その曲線を求める式を用意しておこう。

二点間の最短コースは先ほどの正距方位図法上で直線で表現できる。その直線は終点$(\rho_e,\psi_e)$と媒介変数$t$を用いて、

$$u=t\rho_e\cos\psi_e$$
$$v=t\rho_e\sin\psi_e$$

つまり、

$$\rho=t\rho_e$$
$$\psi=\psi_e$$

である。これは一定方位を維持したまま緯度を一定速度で進行するのは最短コースであるという式である。

これをメルカトル図法上で表すために、一旦球座標に戻そう。
なお先ほど定義した通り$\rho$は$t$に依存した関数なので、ここから先はすべて$t$に依存した関数になることに注意する。

$$x'=\sin\rho\cos\psi$$
$$y'=\sin\rho\sin\psi$$
$$z'=\cos\rho$$

先ほどの回転行列(ユニタリ行列)の逆行列をかけて、元の球座標に戻す。
ユニタリ行列の逆行列は、成分が実数ならば転置と同じなので

$$X=U^TX'$$

これを球座標に戻す。

$$\theta=\arccos z\quad (0\leq\theta\leq\pi)$$
$$\phi=\arccos \frac{x}{\sin\theta+\epsilon}\mathrm{sign}(y)$$

あとは、これをメルカトル図法座標(u,v)に直す

$$u=\phi$$
$$v=\ln\tan\frac{\theta}{2}$$

この手順で得られた関数$u(t),v(t)$がメルカトル図法上での最短距離を表す曲線になる。

これを実際にシミュレーションしてみたいと思うが、それは次回に続く。