人工衛星の軌道やGNSSログなど、XYZの形で表された座標データを、緯度・経度・高度の形に変換する方法について考えてみよう。
一見すると簡単な幾何学の問題に思えるが、地球を「球」ではなく「回転楕円体」として扱おうとすると、途端に話はややこしくなる。
ECEF（Earth Centered Earth Fixed）座標から測地（geodetic）座標への変換は、実は解析的に解くのが難しく、古くから多くの研究者が効率的な解法を競ってきたトピックだ。
今回は、測地座標への変換の基礎と、PROJの実装にも採用されているBowringの手法、そしてVermeilleによる解析的なアプローチについて紹介していきたい。
ECEF座標と測地（geodetic）座標の関係
まず地球表面上の点r s i t e \boldsymbol{r}_\mathrm{site} r site φ r d \varphi_\mathrm{rd} φ rd φ g c \varphi_\mathrm{gc} φ gc
r s i t e = [ R ⊕ cos φ r d cos λ R ⊕ cos φ r d sin λ b ⊕ sin φ r d ] \begin{gather}
\boldsymbol{r}_\mathrm{site} = \left[ \begin{array}{c}
R_\oplus \cos\varphi_\mathrm{rd} \cos\lambda \\
R_\oplus \cos\varphi_\mathrm{rd} \sin\lambda \\
b_\oplus \sin\varphi_\mathrm{rd} \end{array} \right]
\end{gather} r site = R ⊕ cos φ rd cos λ R ⊕ cos φ rd sin λ b ⊕ sin φ rd Figure 1: Spheroidal Earth Geometry and Geodetic Coordinate Parameters.
一方で、geocentric latitude φ g c \varphi_\mathrm{gc} φ gc r s i t e r_\mathrm{site} r site
r s i t e = [ r s i t e cos φ g c cos λ r s i t e cos φ g c sin λ r s i t e sin φ g c ] \begin{gather}
\boldsymbol{r}_\mathrm{site} = \left[ \begin{array}{c}
r_\mathrm{site} \cos\varphi_\mathrm{gc} \cos\lambda \\
r_\mathrm{site} \cos\varphi_\mathrm{gc} \sin\lambda \\
r_\mathrm{site} \sin\varphi_\mathrm{gc} \end{array} \right]
\end{gather} r site = r site cos φ gc cos λ r site cos φ gc sin λ r site sin φ gc これらの式を見比べると、φ r d \varphi_\mathrm{rd} φ rd φ g c \varphi_\mathrm{gc} φ gc
tan φ g c = b ⊕ R ⊕ tan φ r d \begin{gather}
\tan \varphi_\mathrm{gc} = \frac{b_\oplus}{R_\oplus} \tan \varphi_\mathrm{rd}
\end{gather} tan φ gc = R ⊕ b ⊕ tan φ rd 次に、φ r d \varphi_\mathrm{rd} φ rd φ g d \varphi_\mathrm{gd} φ gd r s i t e \boldsymbol{r}_\mathrm{site} r site φ r d \varphi_\mathrm{rd} φ rd λ \lambda λ
d d φ r d r s i t e = [ − R ⊕ sin φ r d cos λ − R ⊕ sin φ r d sin λ b ⊕ cos φ r d ] , d d λ r s i t e = [ − R ⊕ cos φ r d sin λ + R ⊕ cos φ r d cos λ 0 ] . \begin{align}
\frac{d}{d\varphi_\mathrm{rd}} \boldsymbol{r}_\mathrm{site} &= \left[ \begin{array}{c}
-R_\oplus \sin\varphi_\mathrm{rd} \cos\lambda \\
-R_\oplus \sin\varphi_\mathrm{rd} \sin\lambda \\
b_\oplus \cos\varphi_\mathrm{rd} \end{array} \right], \\
\frac{d}{d\lambda} \boldsymbol{r}_\mathrm{site} &=
\left[ \begin{array}{c}
-R_\oplus\cos\varphi_\mathrm{rd} \sin\lambda \\
+R_\oplus\cos\varphi_\mathrm{rd} \cos\lambda \\
0
\end{array} \right]. \\
\end{align} d φ rd d r site d λ d r site = − R ⊕ sin φ rd cos λ − R ⊕ sin φ rd sin λ b ⊕ cos φ rd , = − R ⊕ cos φ rd sin λ + R ⊕ cos φ rd cos λ 0 . で、これをよく見ると、r s i t e \bm{r}_\mathrm{site} r site
n = [ b ⊕ cos φ r d cos λ b ⊕ cos φ r d sin λ R ⊕ sin φ r d ] \begin{gather}
\boldsymbol{n} = \left[ \begin{array}{c}
b_\oplus \cos\varphi_\mathrm{rd} \cos\lambda \\
b_\oplus \cos\varphi_\mathrm{rd} \sin\lambda \\
R_\oplus \sin\varphi_\mathrm{rd}
\end{array} \right]
\end{gather} n = b ⊕ cos φ rd cos λ b ⊕ cos φ rd sin λ R ⊕ sin φ rd 法線方向が分かると、reduced latitude φ r d \varphi_\mathrm{rd} φ rd φ g d \varphi_\mathrm{gd} φ gd
tan φ g d = R ⊕ b ⊕ tan φ r d \begin{gather}
\tan \varphi_\mathrm{gd} = \frac{R_\oplus}{b_\oplus} \tan \varphi_\mathrm{rd}
\end{gather} tan φ gd = b ⊕ R ⊕ tan φ rd 式(3)および式(7)より、geocentric latitude φ g c \varphi_\mathrm{gc} φ gc φ g d \varphi_\mathrm{gd} φ gd
tan φ g c = b ⊕ 2 R ⊕ 2 tan φ g d \begin{gather}
\tan \varphi_\mathrm{gc} = \frac{b_\oplus^2}{R_\oplus^2} \tan \varphi_\mathrm{gd}
\end{gather} tan φ gc = R ⊕ 2 b ⊕ 2 tan φ gd これで、地表面でgeocentricなパラメタからgeodeticなパラメタへ変換することができた。
ただ今知りたいのは、軌道上の点( x , y , z ) (x,y,z) ( x , y , z ) ( x , y , z ) (x,y,z) ( x , y , z ) C ⊕ C_\oplus C ⊕ φ g d \varphi_\mathrm{gd} φ gd
[ x y z ] = [ ( C ⊕ + h e l l p ) cos φ g d cos λ ( C ⊕ + h e l l p ) cos φ g d sin λ ( S ⊕ + h e l l p ) sin φ g d ] \begin{align}
\left[ \begin{array}{c} x \\ y \\ z \end{array} \right] =
\left[ \begin{array}{c}
(C_\oplus + h_\mathrm{ellp}) \cos\varphi_\mathrm{gd} \cos\lambda \\
(C_\oplus + h_\mathrm{ellp}) \cos\varphi_\mathrm{gd} \sin\lambda \\
(S_\oplus + h_\mathrm{ellp}) \sin\varphi_\mathrm{gd}
\end{array} \right]
\end{align} x y z = ( C ⊕ + h ellp ) cos φ gd cos λ ( C ⊕ + h ellp ) cos φ gd sin λ ( S ⊕ + h ellp ) sin φ gd r δ = r s i t e cos φ g c = R ⊕ cos φ r d = C ⊕ cos φ g d r k = r s i t e sin φ g c = R ⊕ 1 − e ⊕ 2 sin φ r d = S ⊕ sin φ g d \begin{align}
r_\delta
&= r_\mathrm{site} \cos \varphi_\mathrm{gc}
= R_\oplus \cos \varphi_\mathrm{rd}
= C_\oplus \cos \varphi_\mathrm{gd} \\
r_\mathrm{k}
&= r_\mathrm{site} \sin \varphi_\mathrm{gc}
% = R_\oplus \frac{b_\oplus}{R_\oplus}\sin \varphi_\mathrm{rd}
= R_\oplus \sqrt{1-e_\oplus^2}\sin \varphi_\mathrm{rd}
= S_\oplus \sin \varphi_\mathrm{gd}
\end{align} r δ r k = r site cos φ gc = R ⊕ cos φ rd = C ⊕ cos φ gd = r site sin φ gc = R ⊕ 1 − e ⊕ 2 sin φ rd = S ⊕ sin φ gd 上式で、cos φ r d \cos \varphi_\mathrm{rd} cos φ rd cos φ r d \cos \varphi_\mathrm{rd} cos φ rd φ g d \varphi_\mathrm{gd} φ gd
tan 2 φ g d = R ⊕ 2 b ⊕ 2 ( 1 cos 2 φ r d − 1 ) 1 cos 2 φ r d = b ⊕ 2 R ⊕ 2 tan 2 φ g d + 1 cos φ r d = 1 b ⊕ 2 R ⊕ 2 tan 2 φ g d + 1 , w h e r e − π 2 ≤ φ r d ≤ π 2 \begin{align}
&\tan^2\varphi_\mathrm{gd} = \frac{R_\oplus^2}{b_\oplus^2} \left(\frac{1}{\cos^2\varphi_\mathrm{rd}} - 1 \right) \\
&\frac{1}{\cos^2\varphi_\mathrm{rd}} = \frac{b_\oplus^2}{R_\oplus^2} \tan^2\varphi_\mathrm{gd} + 1 \\
&\cos\varphi_\mathrm{rd} = \frac{1}{\sqrt{\frac{b_\oplus^2}{R_\oplus^2} \tan^2\varphi_\mathrm{gd} + 1}}, \quad \mathrm{where} \quad -\frac{\pi}{2} \le \varphi_\mathrm{rd} \le \frac{\pi}{2}
\end{align} tan 2 φ gd = b ⊕ 2 R ⊕ 2 ( cos 2 φ rd 1 − 1 ) cos 2 φ rd 1 = R ⊕ 2 b ⊕ 2 tan 2 φ gd + 1 cos φ rd = R ⊕ 2 b ⊕ 2 tan 2 φ gd + 1 1 , where − 2 π ≤ φ rd ≤ 2 π 最終的に、C ⊕ C_\oplus C ⊕
C ⊕ = R ⊕ cos φ g d 1 b ⊕ 2 R ⊕ 2 tan 2 φ g d + 1 = R ⊕ b ⊕ 2 R ⊕ 2 sin 2 φ g d + cos 2 φ g d = R ⊕ ( 1 − e ⊕ 2 ) sin 2 φ g d + cos 2 φ g d = R ⊕ 1 − e ⊕ 2 sin 2 φ g d \begin{align}
C_\oplus &= \frac{R_\oplus}{\cos\varphi_\mathrm{gd}} \frac{1}{\sqrt{\frac{b_\oplus^2}{R_\oplus^2} \tan^2\varphi_\mathrm{gd} + 1}}
= \frac{R_\oplus}{\sqrt{\frac{b_\oplus^2}{R_\oplus^2} \sin^2\varphi_\mathrm{gd} + \cos^2\varphi_\mathrm{gd}}} \notag \\
&= \frac{R_\oplus}{\sqrt{(1-e_\oplus^2) \sin^2\varphi_\mathrm{gd} + \cos^2\varphi_\mathrm{gd}}}
= \frac{R_\oplus}{\sqrt{1 - e_\oplus^2 \sin^2\varphi_\mathrm{gd}}}
\end{align} C ⊕ = cos φ gd R ⊕ R ⊕ 2 b ⊕ 2 tan 2 φ gd + 1 1 = R ⊕ 2 b ⊕ 2 sin 2 φ gd + cos 2 φ gd R ⊕ = ( 1 − e ⊕ 2 ) sin 2 φ gd + cos 2 φ gd R ⊕ = 1 − e ⊕ 2 sin 2 φ gd R ⊕ S ⊕ S_\oplus S ⊕
S ⊕ = R ⊕ ( 1 − e ⊕ 2 ) 1 − e ⊕ 2 sin 2 φ g d \begin{equation}
S_\oplus = \frac{R_\oplus(1-e_\oplus^2)}{\sqrt{1 - e_\oplus^2 \sin^2 \varphi_{\mathrm{gd}}}}
\end{equation} S ⊕ = 1 − e ⊕ 2 sin 2 φ gd R ⊕ ( 1 − e ⊕ 2 ) 数値計算による変換
ここで具体的な変換方法について考えてみよう。
基準となる楕円体としてWGS84を用いることとし、これは長半径（semi-major axis）R ⊕ R_\oplus R ⊕ f ⊕ f_\oplus f ⊕
R ⊕ = 6378.137 k m , f ⊕ = 1 298.257223563 \begin{gather}
R_\oplus = 6378.137~\mathrm{km}, \quad
f_\oplus = \frac{1}{298.257223563}
\end{gather} R ⊕ = 6378.137 km , f ⊕ = 298.257223563 1 その他のパラメタは次のように計算できる。
e ⊕ = 2 f ⊕ − f ⊕ 2 = 0.0818191908426215 , b ⊕ = R ⊕ 1 − e ⊕ 2 = R ⊕ ( 1 − f ⊕ ) = 6356.75231424518 k m . \begin{align}
e_\oplus &= \sqrt{2f_\oplus - f_\oplus^2} = 0.0818191908426215, \\
b_\oplus &= R_\oplus \sqrt{1-e_\oplus^2} = R_\oplus(1-f_\oplus) = 6356.75231424518~\mathrm{km}. \\
\end{align} e ⊕ b ⊕ = 2 f ⊕ − f ⊕ 2 = 0.0818191908426215 , = R ⊕ 1 − e ⊕ 2 = R ⊕ ( 1 − f ⊕ ) = 6356.75231424518 km . 経度に関しては、特に問題なく次のように計算できる。
λ = arctan 2 ( y , x ) \begin{gather}
\lambda = \arctan2 (y,x)
\end{gather} λ = arctan 2 ( y , x ) 緯度と高度に関してはexplicitに解けなさそうなので、数値的に求めることを考えてみる。
まず、仮にφ g d \varphi_\mathrm{gd} φ gd h e l l p h_\mathrm{ellp} h ellp
h e l l p = x 2 + y 2 cos φ g d − C ⊕ = x 2 + y 2 cos φ g d − R ⊕ 1 − e ⊕ 2 sin 2 φ g d h e l l p = z sin φ g d − S ⊕ = z sin φ g d − R ⊕ ( 1 − e ⊕ 2 ) 1 − e ⊕ 2 sin 2 φ g d \begin{align}
h_\mathrm{ellp} &= \frac{\sqrt{x^2+y^2}}{\cos\varphi_\mathrm{gd}} - C_\oplus
= \frac{\sqrt{x^2+y^2}}{\cos\varphi_\mathrm{gd}} - \frac{R_\oplus}{\sqrt{1 - e_\oplus^2 \sin^2\varphi_\mathrm{gd}}} \\
h_\mathrm{ellp} &= \frac{z}{\sin\varphi_\mathrm{gd}} - S_\oplus
= \frac{z}{\sin\varphi_\mathrm{gd}} - \frac{R_\oplus(1-e_\oplus^2)}{\sqrt{1 - e_\oplus^2 \sin^2 \varphi_{\mathrm{gd}}}}
\end{align} h ellp h ellp = cos φ gd x 2 + y 2 − C ⊕ = cos φ gd x 2 + y 2 − 1 − e ⊕ 2 sin 2 φ gd R ⊕ = sin φ gd z − S ⊕ = sin φ gd z − 1 − e ⊕ 2 sin 2 φ gd R ⊕ ( 1 − e ⊕ 2 ) これらが等しくなるように、以下の関係を満たすφ g d \varphi_\mathrm{gd} φ gd
z − x 2 + y 2 tan φ g d + e ⊕ 2 R ⊕ sin φ g d 1 − e ⊕ 2 sin 2 φ g d = 0 \begin{gather}
z - \sqrt{x^2+y^2} \tan \varphi_\mathrm{gd} + \frac{e_\oplus^2 R_\oplus \sin\varphi_\mathrm{gd}}{\sqrt{1 - e_\oplus^2 \sin^2 \varphi_{\mathrm{gd}}}} = 0
\end{gather} z − x 2 + y 2 tan φ gd + 1 − e ⊕ 2 sin 2 φ gd e ⊕ 2 R ⊕ sin φ gd = 0 φ 0 = arctan ( z x 2 + y 2 ) \begin{gather}
\varphi_0 = \arctan \left( \frac{z}{\sqrt{x^2 + y^2}} \right)
\end{gather} φ 0 = arctan ( x 2 + y 2 z ) 初期値はgeocentricなパラメタを用いることとして、
例えばPythonで以下のようなスクリプトを書くと、実際にgeodetic latitudeを求めることができる。
# WGS84 parameters flattening = 1 / 298.257223563 semimajor = 6378.137 eccentricity = np.sqrt( 2 * flattening - flattening ** 2 ) def geodetic_latitude (x, y, z): """ calculate the geodetic latitude # Args: x(ndarray): x coordinate in ECEF, km y(ndarray): y coordinate in ECEF, km z(ndarray): z coordinate in ECEF, km # Returns: latitude(ndarray): geodetic latitude, rad """ params = (x, y, z) latitude = np.arctan(z / np.sqrt(x ** 2 + y ** 2 )) return optimize.fsolve(latitude_equation, latitude, args = params) def latitude_equation (x: np.ndarray, * args) -> np.ndarray: """ equation to be solved # Args: *args(tuple): (x, y, z) # Returns: latitude(ndarray): latitude """ return args[ 2 ] - np.sqrt(args[ 0 ] ** 2 + args[ 1 ] ** 2 ) * np.tan(x) + eccentricity ** 2 * semimajor * np.sin(x) / np.sqrt( 1 - eccentricity ** 2 * np.sin(x) ** 2 ) これでECEF座標からgeodeticなパラメタへの変換の基本は理解できた。
ただこの方法で解を探すと、データ量が増えるにつれて膨大な計算量が必要になってしまうので、より効率的な方法が欲しくなってくる。
Vermeilleの解析的な変換手法
先ほどは、導出した式の複雑さから、解析解を求めることは早々に諦めてしまったが、
実は４次の代数方程式に帰着させることで、Geodetic latitudeと高度を解析的に求めることが出来る。
ここではVermeilleによって提案された、解析的な変換手法について見てみよう[2] 。
まず、次のような変数k k k k > 0 k>0 k > 0
k = Q S P R = h e l l p + C ⊕ − e ⊕ 2 C ⊕ C ⊕ h e l l p = ( k + e ⊕ 2 − 1 ) C ⊕ = k C ⊕ − S ⊕ \begin{align}
&k = \frac{QS}{PR} = \frac{h_\mathrm{ellp} + C_\oplus - e_\oplus^2 C_\oplus}{C_\oplus} \\
&h_\mathrm{ellp} = (k + e_\oplus^2 - 1) C_\oplus = k C_\oplus - S_\oplus
\end{align} k = PR QS = C ⊕ h ellp + C ⊕ − e ⊕ 2 C ⊕ h ellp = ( k + e ⊕ 2 − 1 ) C ⊕ = k C ⊕ − S ⊕ このk k k C ⊕ C_\oplus C ⊕ k k k
sin φ g d = z S ⊕ + h e l l p = z k C ⊕ \begin{align}
\sin\varphi_\mathrm{gd} = \frac{z}{S_\oplus + h_\mathrm{ellp}}
= \frac{z}{k C_\oplus}
\end{align} sin φ gd = S ⊕ + h ellp z = k C ⊕ z C ⊕ 2 = R ⊕ 2 ( 1 − e ⊕ 2 sin 2 φ g d ) = R ⊕ 2 + C ⊕ 2 e ⊕ 2 sin 2 φ g d = R ⊕ 2 + e ⊕ 2 z 2 k 2 \begin{align}
C_\oplus^2 &= \frac{R_\oplus^2}{(1 - e_\oplus^2 \sin^2\varphi_\mathrm{gd})} = R_\oplus^2 + C_\oplus^2 e_\oplus^2 \sin^2\varphi_\mathrm{gd} \notag \\
&= R_\oplus^2 + \frac{e_\oplus^2 z^2}{k^2}
\end{align} C ⊕ 2 = ( 1 − e ⊕ 2 sin 2 φ gd ) R ⊕ 2 = R ⊕ 2 + C ⊕ 2 e ⊕ 2 sin 2 φ gd = R ⊕ 2 + k 2 e ⊕ 2 z 2 これらを用いて、x 2 + y 2 x^2+y^2 x 2 + y 2
x 2 + y 2 = ( h e l l p + C ⊕ ) 2 cos 2 φ g d = ( k + e ⊕ ) 2 C ⊕ 2 ( 1 − sin 2 φ g d ) = ( k + e ⊕ ) 2 C ⊕ 2 ( 1 − z 2 k 2 C ⊕ 2 ) = ( k + e ⊕ ) 2 ( R ⊕ 2 + e ⊕ 2 z 2 k 2 − z 2 k 2 ) \begin{align}
x^2 + y^2
&=(h_\mathrm{ellp} + C_\oplus)^2 \cos^2\varphi_\mathrm{gd}
= (k + e_\oplus)^2 C_\oplus^2 (1 - \sin^2\varphi_\mathrm{gd}) \notag \\
&=(k + e_\oplus)^2 C_\oplus^2 \left( 1 - \frac{z^2}{k^2 C_\oplus^2} \right) \notag \\
&=(k + e_\oplus)^2 \left( R_\oplus^2 + \frac{e_\oplus^2 z^2}{k^2} - \frac{z^2}{k^2}\right)
\end{align} x 2 + y 2 = ( h ellp + C ⊕ ) 2 cos 2 φ gd = ( k + e ⊕ ) 2 C ⊕ 2 ( 1 − sin 2 φ gd ) = ( k + e ⊕ ) 2 C ⊕ 2 ( 1 − k 2 C ⊕ 2 z 2 ) = ( k + e ⊕ ) 2 ( R ⊕ 2 + k 2 e ⊕ 2 z 2 − k 2 z 2 ) x 2 + y 2 ( k + e ⊕ 2 ) 2 + ( 1 − e ⊕ 2 ) z 2 k 2 = R ⊕ 2 \begin{align}
\frac{x^2+y^2}{(k+e_\oplus^2)^2} + \frac{(1 - e_\oplus^2) z^2}{k^2} = R_\oplus^2
\end{align} ( k + e ⊕ 2 ) 2 x 2 + y 2 + k 2 ( 1 − e ⊕ 2 ) z 2 = R ⊕ 2 ここで、p , q p, q p , q
p = x 2 + y 2 R ⊕ 2 , q = 1 − e ⊕ 2 R ⊕ 2 z 2 \begin{align}
p = \frac{x^2 + y^2}{R_\oplus^2}, \quad q = \frac{1-e_\oplus^2}{R_\oplus^2} z^2
\end{align} p = R ⊕ 2 x 2 + y 2 , q = R ⊕ 2 1 − e ⊕ 2 z 2 これらを用いて、k k k
k 4 + 2 e ⊕ 2 k 3 − ( p + q − e ⊕ 4 ) k 2 − 2 e ⊕ 2 q k − e ⊕ 4 q = 0 \begin{align}
k^4 + 2e_\oplus^2 k^3 - (p+q-e_\oplus^4) k^2 - 2e_\oplus^2 q k - e_\oplus^4 q = 0
\end{align} k 4 + 2 e ⊕ 2 k 3 − ( p + q − e ⊕ 4 ) k 2 − 2 e ⊕ 2 q k − e ⊕ 4 q = 0 ここで謎のパラメタu u u u u u
( k 2 + e ⊕ 2 k − u ) 2 − [ ( p + q − 2 u ) k 2 + 2 e ⊕ 2 ( q − u ) k + u 2 + e ⊕ 4 q ] = 0 \begin{align}
(k^2 + e_\oplus^2 k - u)^2 - \left[ (p+q-2u) k^2 + 2e_\oplus^2(q-u)k + u^2 + e_\oplus^4 q \right] = 0
\end{align} ( k 2 + e ⊕ 2 k − u ) 2 − [ ( p + q − 2 u ) k 2 + 2 e ⊕ 2 ( q − u ) k + u 2 + e ⊕ 4 q ] = 0 カギ括弧の中身はk k k ( ⋯ ) 2 (\cdots)^2 ( ⋯ ) 2
e ⊕ 4 ( q − u ) 2 − ( p + q − 2 u ) ( u 2 + e ⊕ 4 q ) = 0 \begin{align}
e_\oplus^4 (q-u)^2 - (p+q-2u)(u^2 + e_\oplus^4 q) = 0
\end{align} e ⊕ 4 ( q − u ) 2 − ( p + q − 2 u ) ( u 2 + e ⊕ 4 q ) = 0 2 u 3 − ( p + q − e ⊕ 4 ) u 2 + e ⊕ 4 p q = 0 \begin{align}
2u^3 - (p + q - e_\oplus^4) u^2 + e_\oplus^4 pq = 0
\end{align} 2 u 3 − ( p + q − e ⊕ 4 ) u 2 + e ⊕ 4 pq = 0 さらに、変数r , s r, s r , s
r = p + q − e ⊕ 4 6 , s = e ⊕ 4 p q 4 r 3 \begin{align}
r = \frac{p+q-e_\oplus^4}{6}, \quad s = e_\oplus^4 \frac{pq}{4r^3}
\end{align} r = 6 p + q − e ⊕ 4 , s = e ⊕ 4 4 r 3 pq ただし、これらの変数は常に正である。
p + q = x 2 + y 2 R ⊕ 2 + ( 1 − e ⊕ 2 ) z 2 R ⊕ 2 = x 2 R ⊕ 2 + y 2 R ⊕ 2 + z 2 b ⊕ 2 > 1 \begin{align}
p + q = \frac{x^2 + y^2}{R_\oplus^2} + \frac{(1-e_\oplus^2) z^2}{R_\oplus^2}
= \frac{x^2}{R_\oplus^2} + \frac{y^2}{R_\oplus^2} + \frac{z^2}{b_\oplus^2} > 1
\end{align} p + q = R ⊕ 2 x 2 + y 2 + R ⊕ 2 ( 1 − e ⊕ 2 ) z 2 = R ⊕ 2 x 2 + R ⊕ 2 y 2 + b ⊕ 2 z 2 > 1 これらを用いて、さらに式(35)の全体を2 r 3 2r^3 2 r 3 u / r u/r u / r
u 3 r 3 − 3 u 2 r 2 − 2 s = 0 \begin{align}
\frac{u^3}{r^3} - 3\frac{u^2}{r^2} - 2s = 0
\end{align} r 3 u 3 − 3 r 2 u 2 − 2 s = 0 ここでu / r u/r u / r t t t t > 0 t>0 t > 0 t = 1 t=1 t = 1 1 + t + 1 t 1+t+\frac{1}{t} 1 + t + t 1 − 2 s -2s − 2 s t t t u r \frac{u}{r} r u 0 < t < 1 0<t<1 0 < t < 1 t > 1 t>1 t > 1 x + 1 x = t x+\frac{1}{x}=t x + x 1 = t
u r = 1 + t + 1 t \begin{align}
\frac{u}{r} = 1 + t + \frac{1}{t}
\end{align} r u = 1 + t + t 1 t 6 − 2 ( 1 + s ) t 3 + 1 = 0 \begin{gather}
t^6 - 2(1+s)t^3 + 1 = 0
\end{gather} t 6 − 2 ( 1 + s ) t 3 + 1 = 0 実際、0 < t < 1 0<t<1 0 < t < 1 t > 1 t>1 t > 1
t 3 = 1 + s ± s ( 2 + s ) t = 1 + s ± s ( 2 + s ) 3 \begin{align}
t^3 &= 1 + s \pm \sqrt{s(2+s)} \\
t &= \sqrt[3]{1 + s \pm \sqrt{s(2+s)}}
\end{align} t 3 t = 1 + s ± s ( 2 + s ) = 3 1 + s ± s ( 2 + s ) で、いずれのt t t u r \frac{u}{r} r u u u u
( k 2 + e ⊕ 2 k − u ) 2 − ( e ⊕ 2 q − u v k + v ) 2 = 0 , w h e r e v = u 2 + e ⊕ 4 q \begin{align}
(k^2 + e_\oplus^2 k - u)^2 - \left( e_\oplus^2 \frac{q-u}{v}k + v \right)^2 = 0, \quad
\mathrm{where}\quad v = \sqrt{u^2 + e_\oplus^4 q}
\end{align} ( k 2 + e ⊕ 2 k − u ) 2 − ( e ⊕ 2 v q − u k + v ) 2 = 0 , where v = u 2 + e ⊕ 4 q ( k 2 + v − u + q v e ⊕ 2 k + v − u ) ( k 2 + v − u + q v e ⊕ 2 k − v − u ) = 0 \begin{align}
\left( k^2 + \frac{v-u+q}{v}e_\oplus^2k + v - u\right)\left( k^2 + \frac{v-u+q}{v}e_\oplus^2k - v - u\right) = 0
\end{align} ( k 2 + v v − u + q e ⊕ 2 k + v − u ) ( k 2 + v v − u + q e ⊕ 2 k − v − u ) = 0 式(44)のひとつ目の括弧内は、v − u , v , q v-u, v, q v − u , v , q k > 0 k>0 k > 0 u + v u+v u + v k > 0 k>0 k > 0
k = u + v + w 2 − w , w h e r e w = e ⊕ 2 u + v − q 2 v \begin{gather}
k = \sqrt{u+v+w^2} - w, \quad \mathrm{where} \quad w = e_\oplus^2\frac{u+v-q}{2v}
\end{gather} k = u + v + w 2 − w , where w = e ⊕ 2 2 v u + v − q これでk k k φ g d \varphi_\mathrm{gd} φ gd h e l l p h_\mathrm{ellp} h ellp
D = k x 2 + y 2 k + e ⊕ 2 , C ⊕ = D 2 + z 2 k \begin{equation}
D = \frac{k\sqrt{x^2+y^2}}{k+e_\oplus^2}, \quad
C_\oplus = \frac{\sqrt{D^2+z^2}}{k}
\end{equation} D = k + e ⊕ 2 k x 2 + y 2 , C ⊕ = k D 2 + z 2 h e l l p = k + e ⊕ 2 − 1 k D 2 + z 2 , φ g d = 2 arctan z D + D 2 + z 2 \begin{equation}
h_\mathrm{ellp} = \frac{k+e_\oplus^2-1}{k} \sqrt{D^2 + z^2}, \quad
\varphi_\mathrm{gd} = 2 \arctan \frac{z}{D+\sqrt{D^2+z^2}}
\end{equation} h ellp = k k + e ⊕ 2 − 1 D 2 + z 2 , φ gd = 2 arctan D + D 2 + z 2 z 元論文[2] には変換に最低限必要な式がリストアップされているので、
それらをもとに簡単なスクリプトを書くと、ECEF座標からgeodeticなパラメタへの変換を実行できる。
# WGS84 parameters flattening = 1 / 298.257223563 semimajor = 6378.137 eccentricity = np.sqrt( 2 * flattening - flattening ** 2 ) def geodetic_vermeille (x, y, z): """ calculate the geodetic latitude and altitude based on H. Vermeille, "Direct transformation from geocentric to geodetic coordinates", 2002, Journal of Geodesy, 76:451-454 # Args: x(ndarray): x coordinate in ECEF, km y(ndarray): y coordinate in ECEF, km z(ndarray): z coordinate in ECEF, km # Returns: lat(ndarray): geodetic latitude, rad lon(ndarray): geodetic longitude, rad h(ndarray): geodetic altitude, km """ p = (x ** 2 + y ** 2 ) / semimajor ** 2 q = ( 1 - eccentricity ** 2 ) * z ** 2 / semimajor ** 2 r = (p + q - eccentricity ** 4 ) / 6 s = eccentricity ** 4 * p * q / ( 4 * r ** 3 ) t = ( 1 + s + np.sqrt(s * ( 2 + s))) ** ( 1 / 3 ) u = r * ( 1 + t + 1 / t) v = np.sqrt(u ** 2 + eccentricity ** 4 * q) w = eccentricity ** 2 * (u + v - q) / ( 2 * v) k = np.sqrt(u + v + w ** 2 ) - w D = k * np.sqrt(x ** 2 + y ** 2 ) / (k + eccentricity ** 2 ) lat = 2 * np.arctan(z / (D + np.sqrt(D ** 2 + z ** 2 ))) lon = np.arctan2(y, x) h = (k + eccentricity ** 2 - 1 ) / k * np.sqrt(D ** 2 + z ** 2 ) return lat, lon, h Bowringの反復的な変換手法
解析解は理論的には美しいものの、コンピュータ上で計算する以上、結局のところ浮動小数点誤差の影響は避けられない。
数値解であったとしても、十分な精度でかつより少ない計算量で解けるならば、実用上はそちらの方が望ましい場合も多い。
実際、PROJ v9.7.0ではBowringによる手法が採用されている[3] 。
再度、(23)を出発点に式変形を行ってみよう。
x 2 + y 2 tan φ g d = z + e ⊕ 2 C ⊕ sin φ g d \begin{equation}
\sqrt{x^2 + y^2} \tan \varphi_\mathrm{gd} = z + e_\oplus^2 C_\oplus \sin \varphi_\mathrm{gd}
\end{equation} x 2 + y 2 tan φ gd = z + e ⊕ 2 C ⊕ sin φ gd (8)および(11)を用いて、φ g d \varphi_\mathrm{gd} φ gd φ r d \varphi_\mathrm{rd} φ rd 1 − e ⊕ 2 = b ⊕ 2 R ⊕ 2 1 - e_\oplus^2 = \frac{b_\oplus^2}{R_\oplus^2} 1 − e ⊕ 2 = R ⊕ 2 b ⊕ 2
x 2 + y 2 R ⊕ b ⊕ tan φ r d = z + e ⊕ 2 C ⊕ R ⊕ S ⊕ 1 − e ⊕ 2 sin φ r d \begin{equation}
\sqrt{x^2 + y^2} \frac{R_\oplus}{b_\oplus} \tan \varphi_\mathrm{rd} = z + e_\oplus^2 C_\oplus \frac{R_\oplus}{S_\oplus} \sqrt{1 - e_\oplus^2} \sin \varphi_\mathrm{rd}
\end{equation} x 2 + y 2 b ⊕ R ⊕ tan φ rd = z + e ⊕ 2 C ⊕ S ⊕ R ⊕ 1 − e ⊕ 2 sin φ rd tan φ r d = z b ⊕ R ⊕ x 2 + y 2 + e ⊕ 2 1 − e ⊕ 2 1 x 2 + y 2 b ⊕ 2 R ⊕ sin φ r d = ( 1 − e ⊕ 2 ) z R ⊕ b ⊕ x 2 + y 2 + e ⊕ 2 R ⊕ x 2 + y 2 sin φ r d \begin{align}
\tan \varphi_\mathrm{rd} &= \frac{z b_\oplus}{R_\oplus\sqrt{x^2 + y^2}} + \frac{e_\oplus^2}{1 - e_\oplus^2} \frac{1}{\sqrt{x^2 + y^2}} \frac{b_\oplus^2}{R_\oplus} \sin \varphi_\mathrm{rd} \notag \\
&= (1 - e_\oplus^2) \frac{z R_\oplus}{b_\oplus\sqrt{x^2 + y^2}} + \frac{e_\oplus^2 R_\oplus}{\sqrt{x^2 + y^2}} \sin \varphi_\mathrm{rd}
\end{align} tan φ rd = R ⊕ x 2 + y 2 z b ⊕ + 1 − e ⊕ 2 e ⊕ 2 x 2 + y 2 1 R ⊕ b ⊕ 2 sin φ rd = ( 1 − e ⊕ 2 ) b ⊕ x 2 + y 2 z R ⊕ + x 2 + y 2 e ⊕ 2 R ⊕ sin φ rd ここで、tan θ = z R ⊕ b ⊕ x 2 + y 2 \tan\theta = \frac{z R_\oplus}{b_\oplus \sqrt{x^2 + y^2}} tan θ = b ⊕ x 2 + y 2 z R ⊕
tan φ r d − ( 1 − e ⊕ 2 ) tan θ − e ⊕ 2 R ⊕ x 2 + y 2 sin φ r d = 0 \begin{equation}
\tan \varphi_\mathrm{rd} - (1 - e_\oplus^2) \tan\theta - \frac{e_\oplus^2 R_\oplus}{\sqrt{x^2 + y^2}} \sin \varphi_\mathrm{rd} = 0
\end{equation} tan φ rd − ( 1 − e ⊕ 2 ) tan θ − x 2 + y 2 e ⊕ 2 R ⊕ sin φ rd = 0 この式を満たすようなtan φ r d \tan \varphi_\mathrm{rd} tan φ rd tan θ \tan\theta tan θ f ( t ) f(t) f ( t ) f ′ ( t ) f'(t) f ′ ( t )
t n + 1 = t n − Δ t , f ( t n ) = f ′ ( t n ) Δ t \begin{equation}
t_{n+1} = t_n - \Delta t, \quad f(t_n) = f'(t_n) \Delta t
\end{equation} t n + 1 = t n − Δ t , f ( t n ) = f ′ ( t n ) Δ t f ( t ) = t − ( 1 − e ⊕ 2 ) tan θ − e ⊕ 2 R ⊕ x 2 + y 2 t 1 + t 2 \begin{equation}
f(t) = t - (1 - e_\oplus^2) \tan\theta - \frac{e_\oplus^2 R_\oplus}{\sqrt{x^2 + y^2}} \frac{t}{\sqrt{1 + t^2}}
\end{equation} f ( t ) = t − ( 1 − e ⊕ 2 ) tan θ − x 2 + y 2 e ⊕ 2 R ⊕ 1 + t 2 t f ′ ( t ) = 1 − e ⊕ 2 R ⊕ x 2 + y 2 ( ( 1 + t 2 ) − 1 2 − t 2 ( 1 + t 2 ) − 3 2 ) = 1 − e ⊕ 2 R ⊕ x 2 + y 2 1 ( 1 + t 2 ) 3 2 \begin{align}
f'(t) &= 1 - \frac{e_\oplus^2 R_\oplus}{\sqrt{x^2 + y^2}} \left( (1 + t^2)^{-\frac{1}{2}} - t^2 (1 + t^2)^{-\frac{3}{2}} \right) \notag \\
&=1 - \frac{e_\oplus^2 R_\oplus}{\sqrt{x^2 + y^2}} \frac{1}{(1 + t^2)^{\frac{3}{2}}}
\end{align} f ′ ( t ) = 1 − x 2 + y 2 e ⊕ 2 R ⊕ ( ( 1 + t 2 ) − 2 1 − t 2 ( 1 + t 2 ) − 2 3 ) = 1 − x 2 + y 2 e ⊕ 2 R ⊕ ( 1 + t 2 ) 2 3 1 これより具体的な推定式は次のように表される。
tan φ r d ≈ tan θ − f ( tan θ ) f ′ ( tan θ ) = tan θ − e ⊕ 2 tan θ − e ⊕ 2 R ⊕ x 2 + y 2 sin θ 1 − e ⊕ 2 R ⊕ x 2 + y 2 cos 3 θ = tan θ − e ⊕ 2 tan θ x 2 + y 2 − e ⊕ 2 R ⊕ sin θ x 2 + y 2 − e ⊕ 2 R ⊕ cos 3 θ = ( 1 − e ⊕ 2 ) tan θ x 2 + y 2 − e ⊕ 2 R ⊕ cos 2 θ sin θ + e ⊕ 2 R ⊕ sin θ x 2 + y 2 − e ⊕ 2 R ⊕ cos 3 θ = b ⊕ R ⊕ z + e ⊕ 2 R ⊕ sin 3 θ x 2 + y 2 − e ⊕ 2 R ⊕ cos 3 θ \begin{align}
\tan \varphi_\mathrm{rd} &\approx \tan\theta - \frac{f(\tan\theta)}{f'(\tan\theta)} \notag \\
&= \tan\theta - \frac{e_\oplus^2 \tan\theta - \frac{e_\oplus^2 R_\oplus}{\sqrt{x^2 + y^2}} \sin\theta}{1 - \frac{e_\oplus^2 R_\oplus}{\sqrt{x^2 + y^2}} \cos^3\theta} \notag \\
&= \tan\theta - \frac{e_\oplus^2 \tan\theta \sqrt{x^2 + y^2} - e_\oplus^2 R_\oplus \sin\theta}{\sqrt{x^2 + y^2} - e_\oplus^2 R_\oplus \cos^3\theta} \notag \\
&= \frac{(1 - e_\oplus^2) \tan\theta \sqrt{x^2 + y^2} - e_\oplus^2 R_\oplus \cos^2\theta\sin\theta + e_\oplus^2 R_\oplus \sin\theta}{\sqrt{x^2 + y^2} - e_\oplus^2 R_\oplus \cos^3\theta} \notag \\
&= \frac{\frac{b_\oplus}{R_\oplus}z + e_\oplus^2 R_\oplus \sin^3\theta}{\sqrt{x^2 + y^2} - e_\oplus^2 R_\oplus \cos^3\theta}
\end{align} tan φ rd ≈ tan θ − f ′ ( tan θ ) f ( tan θ ) = tan θ − 1 − x 2 + y 2 e ⊕ 2 R ⊕ cos 3 θ e ⊕ 2 tan θ − x 2 + y 2 e ⊕ 2 R ⊕ sin θ = tan θ − x 2 + y 2 − e ⊕ 2 R ⊕ cos 3 θ e ⊕ 2 tan θ x 2 + y 2 − e ⊕ 2 R ⊕ sin θ = x 2 + y 2 − e ⊕ 2 R ⊕ cos 3 θ ( 1 − e ⊕ 2 ) tan θ x 2 + y 2 − e ⊕ 2 R ⊕ cos 2 θ sin θ + e ⊕ 2 R ⊕ sin θ = x 2 + y 2 − e ⊕ 2 R ⊕ cos 3 θ R ⊕ b ⊕ z + e ⊕ 2 R ⊕ sin 3 θ 求めたreduced latitude φ r d \varphi_\mathrm{rd} φ rd φ g d \varphi_\mathrm{gd} φ gd
φ g d = arctan ( R ⊕ b ⊕ tan φ r d ) \begin{gather}
\varphi_\mathrm{gd} = \arctan \left(\frac{R_\oplus}{b_\oplus} \tan \varphi_\mathrm{rd}\right)
\end{gather} φ gd = arctan ( b ⊕ R ⊕ tan φ rd ) Reference
David A. Vallado, Fundamentals of Astrodynamics and Applications Fourth Edition, 2013, Microsoft Press
H. Vermeille, “Direct transformation from geocentric to geodetic coordinates”, 2002, Journal of Geodesy, 76:451-454, doi: 10.1007/s00190-002-0273-6 .
Bowring, B. R. (1976). TRANSFORMATION FROM SPATIAL TO GEOGRAPHICAL COORDINATES. Survey Review, 23(181), 323–327. doi: 10.1179/sre.1976.23.181.323
Figure 1: Spheroidal Earth Geometry and Geodetic Coordinate Parameters.