Edited at

JGD2011の座標系にtowgs84が無いとかそもそもtowgs84って何やねん


はじめに

MySQL 8.0.13 の ST_Transform()を試すで「有識者からの情報に期待したいと思います」とのことですので、知識はそんなにないのですが(意識はある(混濁していない))、コメントでもと思ったけど、記事にするとふやかせるので、記事にします。

ですから、これを独立した文書と考えずに、sakaikさんの記事と併せて読んでいただくと趣深くなるかと思います。


towgs84の各要素の意味

towgs84は、WGS84以外の測地系からWGS84に変換するためのパラメータです。3次元直交空間を想定して、平行移動、拡大縮小、回転を行い、Tokyoデータム等からWGS84に無理やり当てはめます。

なお、JGD2000やJGD2011は、WGS84と結構近いので、特に問題がないなら(アンダーメートルの精度が要るほどじゃないなら)、towgs84=[0,0,0,0,0,0,0]として、変換しません。


参考文献

飛田幹男(2002)「世界測地系と座標変換」,日本測量協会,東京,174pp.


それぞれの意味を数式で示すと

towgs84=[t_x,t_y,t_z,r_x,r_y,r_z,d]

tは平行移動、dは倍率、rは回転(ラジアン)です。

回転rの添え字(x, y, z)は、それぞれ、どの軸で回転させるかを指定しています。x軸で回転させると、y軸とz軸からなる平面がぐるっと回転し、y軸で回転させると、z軸とx軸、z軸ならx軸とy軸、といったふういなります。

\left(

\begin{array}{c}
x' \\
y' \\
z'
\end{array}
\right)
=
\left(
\begin{array}{c}
t_1 \\
t_2 \\
t_3
\end{array}
\right)
+
D
\left(
\begin{array}{c}
x \\
y \\
z
\end{array}
\right)
+
R_x R_y R_z \left(
\begin{array}{c}
x \\
y \\
z
\end{array}
\right)

とりあえず、回転行列は出していませんが、まあ計算すると面倒くさいのでお勧めしません。出典元では、r(ラジアン)が非常に小さい数として、微分値で近似したうえ、各要素の2次項以降を捨てています。


訂正 2019-01-07

以前、towgs84=[t_x,t_y,t_z,d,r_x,r_y,r_z]と書いていましたが、座標系変換ライブラリProj.4では、towgs84=[t_x,t_y,t_z,r_x,r_y,r_z,d]となっていました。前者は、Proj.4では誤りです。

https://github.com/OSGeo/proj.4/blob/d0506e19a71888f7f0c3aa8618d919624e754c4d/src/pj_transform.c によると

#define Dx_BF (defn->datum_params[0])

#define Dy_BF (defn->datum_params[1])
#define Dz_BF (defn->datum_params[2])
#define Rx_BF (defn->datum_params[3])
#define Ry_BF (defn->datum_params[4])
#define Rz_BF (defn->datum_params[5])
#define M_BF (defn->datum_params[6])
...
int pj_geocentric_to_wgs84( PJ *defn,
...
{
...
x_out = M_BF*( x[io] - Rz_BF*y[io] + Ry_BF*z[io]) + Dx_BF;

これらよりdに該当するのは、Proj.4ではM_BFで、マクロで配列(0始まり)の6番目のパラメータです。

以上より、以前の7パラメータの順序は間違いでした。お詫びして訂正します。


実質3パラメータの場合

実質3パラメータの場合、つまり次のような場合を考えます。

towgs84=[t_x,t_y,t_z,0,0,0,0]

\begin{align}

\left(
\begin{array}{c}
x' \\
y' \\
z'
\end{array}
\right)
&=
\left(
\begin{array}{c}
t_x \\
t_y \\
t_z
\end{array}
\right)
+
0
\left(
\begin{array}{c}
x \\
y \\
z
\end{array}
\right)
+
E E E \left(
\begin{array}{c}
x \\
y \\
z
\end{array}
\right) \\
&=
\left(
\begin{array}{c}
t_x \\
t_y \\
t_z
\end{array}
\right)
+
\left(
\begin{array}{c}
x \\
y \\
z
\end{array}
\right) \\

\end{align}

ちゃんとシフトしました。


全部0で無変換になるか

ふと、ちゃんと全部0ならx,y,zの値が全て変換されないようになっているか、見てみたくなったので、見てみました。

towgs84=[0,0,0,0,0,0,0]

とすると、

\begin{align}

\left(
\begin{array}{c}
x' \\
y' \\
z'
\end{array}
\right)
&=
\left(
\begin{array}{c}
0 \\
0 \\
0
\end{array}
\right)
+
0
\left(
\begin{array}{c}
x \\
y \\
z
\end{array}
\right)
+
E E E \left(
\begin{array}{c}
x \\
y \\
z
\end{array}
\right) \\
&=
\left(
\begin{array}{c}
x \\
y \\
z
\end{array}
\right) \\

\end{align}

大丈夫でしたね。


なぜ JGD2011には TOWGS84 の記述がないのか

towgs84があると、座標系変換で測地系(というか基準回転楕円対面)を変更する場合に、WGS84を媒介にして、相互変換が可能となります。

EPSG:6668にはtowgs84がなく、MySQLではエラーが出ます。

なんでtowgs84が無いのでしょうか?とりあえず[0,0,0,0,0,0,0]としておけばいいのに?と思うこととでしょう。


知りません。わかりません。

私も知りません。わかりません。

少なくともProj.4epsgファイルでも(それを基にしているPostGISのspatial_ref_sysビューでも)EPSG:6668には、towgs84がありませんでした。

原典と言える https://www.epsg-registry.org/ では、そもそも全てのエントリにtowgs84がありません。

恐いのが、上書きされるかもしれない点です。たぶん、どこからかデータベースを引っ張ってきてMySQLに反映させてるのだと思います(独自管理は相当つらいです)。で、他のデータベースが更新されると、それにあわせてMySQLにも反映させるでしょうけど、この際にtowgs84を追加したのを消される可能性があります。私のようなおっちょこちょいでないならいいだけかも知れないのですが、ちょっと危険かな、という気持ちです。


Proj.4系ではどうなるかというと

db=# SELECT ST_AsEWKT(ST_Transform('SRID=4301;POINT(135 35)'::GEOMETRY, 6668));

st_asewkt
-------------------------
SRID=6668;POINT(135 35)
(1 )

Tokyo地理座標系からJGD2011地理座標系への変換が、スルーされて、かつエラーも警告もありません。ならむしろエラーで撥ねてもらった方がありがたいやないの、と思います。


おわりに

MySQL 8.0.13 の ST_Transform()を試すに対するコメントめいた投稿となりました。

確かFOSS4G Tokyoの懇親会でノートパソコン広げてyyamasaki1さんを縛り付けるという悪事を働いた際に「JGD2011にtowgs84が無い問題」にたまたま接触して、でもProj.4とかの挙動があったので、バグレポートは止めてください、と言ってた記憶があります。

PostGISユーザたる私としては、Tokyo(WGS84に近い測地系でない)の座標系がからんでくる場合に限って、JGD2011系は使わずJGD2000に変換したあとST_SetSRID()で済むから別にいいやとか思ってたりしてましたが、放置はマズいよなー、と。実際は放置しますけどね。


本記事のライセンス

クリエイティブ・コモンズ・ライセンス

この記事は クリエイティブ・コモンズ 表示 4.0 国際 ライセンスの下に提供されています。