JGD2011からのST_Transform()、JGD2011へのST_Transform()

この日記は RDBMS-GIS Advent Calendar20日目の枠です。


 一昨日に書いた MySQL 8.0.13 の ST_Transform()を試す - sakaikの日々雑感~(T)編 に対する返歌を 有意識者の boiledorange73 さんにいただいた(JGD2011の座標系にtowgs84が無いとかそもそもtowgs84って何やねん - Qiita)ので、それに対する恋文返しです。実際には、boiledorange73 さんへのお返事になっているというよりも、boiledorange73 さんにいただいた恋文をこっそりアレンジして他の人に送りつけているような感じと言ったほうが近いかもしれません。

TOWGS84 の 7つの数字の意味

 前の日記で、

(7個ある数字の意味については今後調べていきたいと思います)

 と書きましたが、 boiledorange73 さんが解説してくださいました。ありがとうございます。

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

 で、tが各方向への並行移動量、d が縮尺に関する値、rが回転に関する値であるということです。

SRID:4301(Tokyo)で 、 TOWGS84[-147,506,687,0,0,0,0] と定義されていたのは、最初の3つのみに値が指定されていて残りの4つはゼロであることから、回転や縮尺の変更は伴わずに単純な並行移動のみで変換ができるということをあらわしています。

x,y,z 方向への平行移動の値は、こんなものか?の確認

 この定義によると、x, y, z 方向への並行移動はそれぞれ、-147, 506, 687 とされています。同じ緯度経度の数字とした場合、その指し示す位置は、Tokyo測地系からJGD2000 に向かっておおよそ東南東の方向に約600mほどズレることは、既に知識として持っているので*1、感覚的にこれに合致するかを確認することにします。
 x, y, z と示されていますが、x=-147, y=687 であることから、xが南北方向(北方向をプラス)、yが東西方向(東方向をプラス) に関わる値であると推測できます。
 では、z の 687 というのはどうでしょうか。これも、Tokyo測地系で使用しているベッセル楕円体と、JGD2000が採用しているGRS80 とを比べることで分かります。両者の極半径の差は 674m、赤道半径の差は 740m ですので、まぁなんとなく、論外な値でもなさそうだなと見えます。それっぽいので、終了としたいと思います。(緯度経度の話の中には z軸は出てこない)

行列はキライ

 私は行列は嫌いです。ラーメンは大好きですが、行列してまで食べません。
ということで参照先の紹介だけ。
boiledorange73さんの式の中で出てきた 回転の計算式の R_x, R_y, R_z はそれぞれ、x, y, z方向に関係する 3x3 の行列です。これについては、boiledorange73さんも紹介していた書籍
世界測地系と座標変換
の第2章に詳しく出ているので、興味のある方は参考になさってください。この本、非常に軽妙に、知りたい情報を書いてくれているのですが、私のようなこの分野の素人にとっては、もう少し丁寧に説明して欲しいところもあるというのが率直な感想ですが、3ヶ月前に意味もわからなかった部分が後に少し理解が進んでいることが把握できるなど、自分の成長を感じられる部分もあるのが面白いです。

拡大縮小のパラメタは?

 boiledorange73さんの示してくれた変換式をよく見ると、拡大縮小のパラメタ D を用いる項が目に付きます。「倍率」と示されていますが、倍率であるなら、変化なしの場合に 1 であるべきで、今回のような 0 という設定値には違和感があります。調べてみました。
 正確には、このパラメタ D は、「スケールファクタ(スケール補正量)」と言って、補正なしならゼロ、補正ありなら プラスいくつ、またはマイナスいくつ という値を指定するようです。なお、値の単位は  10^{-8} です。
 

結局 JGD2011 からの/への 座標変換はどうしたらいいの

 JGD2011 と JGD2000 は、地理座標系に於いては同一(=7つの変換パラメタがすべてゼロ)と考えて良いので、boiledorange73さんのエントリでも最後にちらっと書かれていたように、JGD2000を経由して(JGD2011なんだけどJGD2000の値だと見なして)変換処理を行っても問題なさそうです。

一昨日の私のエントリのデータを使って試してみます。

mysql> select ST_AsText(g), ST_SRID(g), ST_AsText( ST_Transform(g, 4301) ) FROM g1  WHERE ST_SRID(g)=6668;
ERROR 3743 (22S00): Transformation from SRID 6668 is not supported. The spatial reference system has no TOWGS84 clause.

 この行のgの値は (towgs84がない)SRID:6668(JGD2011)であるためにエラーになったので、 (変換ではなくそのままの値でSRIDだけ変更する関数である)ST_SRID() を g に噛ますことで、JGD2000ということにして、変換を実施できます。

mysql> select ST_AsText(g), ST_SRID(g), ST_AsText( ST_Transform( ST_SRID(g, 4612) , 4301) ) FROM g1  WHERE ST_SRID(g)=6668;
+---------------+------------+----------------------------------------------------+
| ST_AsText(g)  | ST_SRID(g) | ST_AsText( ST_Transform( ST_SRID(g,4612) , 4301) ) |
+---------------+------------+----------------------------------------------------+
| POINT(35 135) |       6668 | POINT(34.996751668769726 135.00278101474854)       |
+---------------+------------+----------------------------------------------------+


 逆に、Tokyo測地系で与えられた緯度経度を、JGD2011 へと変換したい場合の例は、以下。そのまま 4301(Tokyo)を 6668(JGD2011)へ変換するとエラーになるので、一旦 4612(JGD2000)へと変換し、結果値を 6668(JGD2011)と見なす指示を行うことで、TokyoからJGD2011への変換を実現可能です。

mysql> select ST_AsText(g), ST_SRID(g), ST_AsText( ST_Transform(g, 6668) ) FROM g1  WHERE ST_SRID(g) = 4301;
ERROR 3744 (22S00): Transformation to SRID 6668 is not supported. The spatial reference system has no TOWGS84 clause.

mysql> select ST_AsText(g), ST_SRID(g), ST_AsText(ST_SRID( ST_Transform(g, 4612), 6668) ) FROM g1  WHERE ST_SRID(g) = 4301;
+---------------+------------+---------------------------------------------------+
| ST_AsText(g)  | ST_SRID(g) | ST_AsText(ST_SRID( ST_Transform(g, 4612), 6668) ) |
+---------------+------------+---------------------------------------------------+
| POINT(35 135) |       4301 | POINT(35.003247866817325 134.99721914450956)      |
+---------------+------------+---------------------------------------------------+

 

ST_SRID() 関数って?

 ST_SRID() の2種類の使い方を、何の説明のなく出したので、混乱した人もいるかもしれません。
ST_SRID() 関数は、引数がひとつ(GEOMETRY型)の場合には、それがあらわされているSRIDの数字を返し、第2引数に数字(SRID)を与えた時には、そのSRIDへと変換した GEOMETRY を返す、という挙動をします。
 単に情報を取得したい場合と、変換したい場合とで同じ関数を使うので、慣れるまではちょっとややこしいと私は感じました。

MySQL :: MySQL 8.0 Reference Manual :: 12.16.7.1 General Geometry Property Functions

ということで

 アドヴェントカレンダー用エントリにする程でもなかったかもしれませんが、本日が終わりそうなタイミング(というかこれを書いているまさに今、本日が終わって新しい本日になってしまった)で20日目の枠が空いたままだったので、アドヴェント参加作品(笑)とさせていただきました。
 
f:id:sakaik:20181220121034j:plain
写真は地球。奥のほうにあるのが月。右奥は、日本近郊を切り出した曲面。日本列島が球面の上に乗っていることを体感できます@地図と測量の科学館@つくば

*1:追記:600mではなくて400mでした。指摘くださった ozo360さん、ありがとうございます。600mも差が出る理由はモデルとしている楕円体の違い(ベッセル楕円体とGRS80)で変換するからかなと想像しましたが、今の私にはまだ理解できず、そのうち出直してきます