斜軸変換についてのメモ

GIS の座標変換ライブラリである Proj には、汎用斜軸変換を行うことができる `+proj=ob_tran` が用意されている。

## はじめに

まず単位球面（$a = b = 1, e = 0$）において、緯度経度座標 $(\varphi, \lambda)$ と直交座標 $(x, y, z)$ とのあいだには以下の関係が成り立つ。

```math
\begin{pmatrix}
x \\
y \\
z
\end{pmatrix} =
\begin{pmatrix}
\cos \varphi \cos \lambda \\
\cos \varphi \sin \lambda \\
\sin \varphi
\end{pmatrix}
```

## 球の回転

いま $xyz$ 座標系に存在する単位球面上の任意の点 $P$ について、球体ごと $z\text{-}y\text{-}z$ の順序で回転を行った結果を考える。

すなわち $z$ 軸を軸として $\theta_1$ だけ回転させ、 $y'$ 軸を軸として $\theta_2$ だけ回転させ、 $z''$ 軸を軸として $\theta_3$ だけ回転させる。そうして得られた $x'''y'''z'''$ 座標系を考える。

点 $P$ の、回転する前の $xyz$ 座標系における座標値 $\boldsymbol{P}_s$ を次のようにおく。

```math
\boldsymbol{P}_s = \begin{pmatrix}
x_s \\
y_s \\
z_s
\end{pmatrix}
```

点 $P$ は座標軸とともに一緒に回転するため、回転したあとの点 $P'''$ の $x'''y'''z'''$ 座標系における座標値 $\boldsymbol{P}'''$ は $\boldsymbol{P}_s$ と等しい。ただし、 $x'''y'''z'''$ 座標系の基底は
$\begin{pmatrix}\boldsymbol{\mathrm{e}}'''_x \  \boldsymbol{\mathrm{e}}'''_y \  \boldsymbol{\mathrm{e}}'''_z \end{pmatrix}$ であり、点 $P'''$ の位置ベクトルは $\begin{pmatrix}\boldsymbol{\mathrm{e}}'''_x \  \boldsymbol{\mathrm{e}}'''_y \  \boldsymbol{\mathrm{e}}'''_z \end{pmatrix} \boldsymbol{P}'''$ である。

$x'''y'''z'''$ 座標系の基底は、 $xyz$ 座標系の基底を回転させたものであるため、次のように計算され、回転行列 $\boldsymbol{R}$ とおく。

```math
\begin{array}{ll}
\begin{pmatrix}
\boldsymbol{\mathrm{e}}'''_x \ 
\boldsymbol{\mathrm{e}}'''_y \ 
\boldsymbol{\mathrm{e}}'''_z
\end{pmatrix} & \hspace{-0.7em} =
\begin{pmatrix}
\boldsymbol{\mathrm{e}}_x \ 
\boldsymbol{\mathrm{e}}_y \ 
\boldsymbol{\mathrm{e}}_z
\end{pmatrix} \boldsymbol{R}_z (\theta_1) \boldsymbol{R}_y (\theta_2) \boldsymbol{R}_z (\theta_3) \\
& \hspace{-0.7em} = \begin{pmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{pmatrix}\begin{pmatrix}
c_1 & -s_1 & 0 \\
s_1 & c_1 & 0 \\
0 & 0 & 1
\end{pmatrix}\begin{pmatrix}
c_2 & 0 & s_2 \\
0 & 1 & 0 \\
-s_2 & 0 & c_2
\end{pmatrix}\begin{pmatrix}
c_3 & -s_3 & 0 \\
s_3 & c_3 & 0 \\
0 & 0 & 1
\end{pmatrix} \\
& \hspace{-0.7em} = \begin{pmatrix}
c_1 c_2 c_3 - s_1 s_3 & -c_1 c_2 s_3 - s_1 c_3 & c_1 s_2 \\
s_1 c_2 c_3 + c_1 s_3 & -s_1 c_2 s_3 + c_1 c_3 & s_1 s_2 \\
-s_2 c_3 & s_2 s_3 & c_2
\end{pmatrix} \\
& \hspace{-0.7em} = \boldsymbol{R}
\end{array}
```

次に、点 $P$ を回転した結果である点 $P'''$ に対し、 $xyz$ 座標系における座標値 $\boldsymbol{P}_d$ を考える。

点 $P'''$ は $xyz$ 座標系で表現しても、 $x'''y'''z'''$ 座標系で表現しても、空間的な位置は当然同一である。このため、それぞれの座標系の基底を用い

```math
\begin{pmatrix}\boldsymbol{\mathrm{e}}'''_x \  \boldsymbol{\mathrm{e}}'''_y \  \boldsymbol{\mathrm{e}}'''_z \end{pmatrix} \boldsymbol{P}'''　= \begin{pmatrix}\boldsymbol{\mathrm{e}}'''_x \  \boldsymbol{\mathrm{e}}'''_y \  \boldsymbol{\mathrm{e}}'''_z \end{pmatrix} \boldsymbol{P}_s = \begin{pmatrix}\boldsymbol{\mathrm{e}}_x \  \boldsymbol{\mathrm{e}}_y \  \boldsymbol{\mathrm{e}}_z \end{pmatrix} \boldsymbol{P}_d
```

すなわち

```math
\boldsymbol{P}_d = \boldsymbol{R} \boldsymbol{P}_s \\
\begin{pmatrix}
x_d \\
y_d \\
z_d
\end{pmatrix} = \begin{pmatrix}
\cos \varphi_d \cos \lambda_d \\
\cos \varphi_d \sin \lambda_d \\
\sin \varphi_d
\end{pmatrix} = \boldsymbol{R} \begin{pmatrix}
x_s \\
y_s \\
z_s
\end{pmatrix} = \boldsymbol{R} \begin{pmatrix}
\cos \varphi_s \cos \lambda_s \\
\cos \varphi_s \sin \lambda_s \\
\sin \varphi_s
\end{pmatrix}
```

となる。

:::note warn
対象地点の空間的位置は変化せず、座標系（座標軸）だけが回転するのではなく、座標系とともに物体（地球）が回転を行う。このため**元の座標系**は $x'''y'''z'''$ 座標系となる。

一方、回転により点 $P$ の空間的な位置は変化し、点 $P'''$ に移動する。この回転が行われたあとの $xyz$ 座標系を**斜軸変換後の地球**と考え、投影を行う。
:::

## 北極の回転

元の座標系における北極の座標値 $\boldsymbol{P}_{Ns}$

```math
\boldsymbol{P}_{Ns} = \begin{pmatrix}
0 \\
0 \\
1
\end{pmatrix}
```

は回転により

```math
\boldsymbol{P}_{Nd} = \boldsymbol{R} \boldsymbol{P}_{Ns} = \begin{pmatrix}
\cos \theta_1 \sin \theta_2 \\
\sin \theta_1 \sin \theta_2 \\
\cos \theta_2
\end{pmatrix}
```

となる。ここで $\frac{\pi}{2} - \theta_2 = \varphi_N$ 、 $\theta_1 = \lambda_N$ と置くと

```math
\boldsymbol{P}_{Nd} = \boldsymbol{R} \boldsymbol{P}_{Ns} = \begin{pmatrix}
\cos \varphi_N \cos \lambda_N \\
\cos \varphi_N \sin \lambda_N \\
\sin \varphi_N
\end{pmatrix}
```

となり、回転後の点 $P'''_N$ が $xyz$ 座標系、すなわち**回転後の座標系**において、緯度 $\varphi_N$ 、経度 $\lambda_N$ の位置に相当する地点であることがわかる。なお $\theta_2$ は余緯度に相当する。

Proj の `ob_tran` のパラメータである `o_lat_p` と `o_lon_p` はこの値を指す。

## 式の変形

```math
\begin{array}{rcrl}
z_d & \hspace{-0.7em} = & \hspace{-0.7em} \sin \varphi_d & \hspace{-0.7em} = \sin \varphi_N \sin \varphi_s + \cos \varphi_N \cos \varphi_s \cos \left( \lambda_s + \theta_3 \right) \\
- \sin \lambda_N x_d + \cos \lambda_N y_d & \hspace{-0.7em} = & \hspace{-0.7em} \cos \varphi_d \sin \left( \lambda_d - \lambda_N \right) & \hspace{-0.7em} = \cos \varphi_s \sin \left( \lambda_s + \theta_3 \right)
\end{array}
```

ここで $-\theta_3 = \lambda_0$ とおくと、

```math
\begin{array}{rl}
\sin \varphi_d & \hspace{-0.7em} = \sin \varphi_N \sin \varphi_s + \cos \varphi_N \cos \varphi_s \cos \left( \lambda_s - \lambda_0 \right) \\
\sin \left( \lambda_d - \lambda_N \right) & \hspace{-0.7em} = \cos \varphi_s \sin \left( \lambda_s - \lambda_0 \right) \  / \  \cos \varphi_d
\end{array}
```

となり、記号の差はあるものの [Map Projections -- a Working Manual](https://pubs.er.usgs.gov/publication/pp1395) の式 (5-7) と (5-8) に一致する。

## あとがき

座標系の変換だもんで軸のみを回転させてみたり、 $xyz$ 座標系が変換元で $x'''y'''z'''$ 座標系が変換後と思い込んでて、 `o_lat_p` と `o_lon_p` の扱いがよくわからなくなったり、いろいろ混乱したので。

この記事に書いた内容も本当にあってるのか、よくわからん。