この日記は RDBMS-GIS Advent Calendar の 18日目の枠です。
MySQLの残念な ST_Transform()
ST_Transform() という関数があります。測地系を変換できる機能です。Tokyo測地系(別名日本測地系)で記述されている緯度経度を、JGD2000/2011 に変換したり、あるいは地理座標系の JGD2011の緯度経度を JGD2011の平面直角座標系の9系の座標*1に変換できたりするもの、だと思います。
思います、と書いたのは、MySQLでは後者について、バージョン 8.0.13 の時点ではまだ対応してないからです。そもそも、この ST_Transform()関数、8.0.13が出る前から MySQLのリファレンスマニュアルには記載されていて、でも実際には実装すらされていなかったので私は「蕎麦屋の出前詐欺*2」と呼んでいたものです。このたび限定的な機能ながら、 8.0.13 で実装されましたが、地理座標系どうしの変換の一部にのみ対応しているのが、現在の状況です。
本エントリでは、MySQL 8.0.13 で、ST_Transform() 関数を試してみた結果を紹介します。
データの用意
様々な地理座標系であらわされた「北緯35度、東経135度」のデータを登録しておきます。
CREATE TABLE g1 (g GEOMETRY);
INSERT INTO g1 VALUES ( ST_GeomFromText('POINT(35 135)', 6668));
INSERT INTO g1 VALUES ( ST_GeomFromText('POINT(35 135)', 4612));
INSERT INTO g1 VALUES ( ST_GeomFromText('POINT(35 135)', 4326));
INSERT INTO g1 VALUES ( ST_GeomFromText('POINT(35 135)', 4301));
mysql> SELECT ST_AsText(g), ST_SRID(g) FROM g1;
+---------------+------------+
| ST_AsText(g) | ST_SRID(g) |
+---------------+------------+
| POINT(35 135) | 6668 | ← JGD2011
| POINT(35 135) | 4612 | ← JGD2000
| POINT(35 135) | 4326 | ← WGS84
| POINT(35 135) | 4301 | ← Tokyo
+---------------+------------+いうまでもないことですが、これらの「35度135度」の点は、兵庫県内の、それぞれ異なる場所を指し示しています（実際には Tokyo 測地系のみが大きく離れていて、あとの３つはほぼ同じ場所を指している）
とりあえず実行してみる
登録した４件のデータに対して、ST_Transform() 関数を噛ませてみましょう。
まず、SRID:6668（JGD2011）への変換を試みます。
mysql> SELECT ST_AsText( ST_Transform(g, 6668) ), ST_SRID(g) FROM g1; ERROR 3744 (22S00): Transformation to SRID 6668 is not supported. The spatial reference system has no TOWGS84 clause.
あれれ。「6668への変換はサポートされてないよ」だそうです。
うむ。では、SRID:4326（WGS84）への変換を試みてみましょう。
mysql> SELECT ST_AsText( ST_Transform(g, 4326) ), ST_SRID(g) FROM g1; ERROR 3743 (22S00): Transformation from SRID 6668 is not supported. The spatial reference system has no TOWGS84 clause.
少しエラーコードとエラーメッセージの内容が変わりました。「6668からの変換はサポートしてないよ」だそうです。
ということで、6668はダメらしいと分かったので、6668以外のものを 4326への変換を試みることにします。残念な6668。
mysql> SELECT ST_AsText( ST_Transform(g, 4326) ), ST_SRID(g) org_srid FROM g1 WHERE ST_SRID(g)<>6668; +---------------------------------------------+----------+ | ST_AsText( ST_Transform(g, 4326) ) | org_srid | +---------------------------------------------+----------+ | POINT(35 135) | 4612 | | POINT(35 135) | 4326 | | POINT(35.00324786593045 134.99721914450956) | 4301 | +---------------------------------------------+----------+
無事実行できて、結果が返ってきました。JGD2000からWGS84への変換は値に変化がありませんが、TokyoからWGS84への変換は、何やら変化があった感じです。この結果は、『Tokyo測地系が「東経135度」と思っている場所は、WGS84に言わせると「東経134.997度」』ということを表しています。（経度だけに着目して語ると）Tokyo測地系の東経135度は明石市役所付近を通り、WGS84やJGD20xxの135度はその東側約400m（天文台の130mくらい西）を通っているという知識を適用すると、これは感覚的には納得できる結果です。
特に意味はないけど、調子に乗って SRID:4269 にも変換してみちゃう。
mysql> SELECT ST_AsText( ST_Transform(g, 4269) ), ST_SRID(g) org_srid FROM g1 WHERE ST_SRID(g)<>6668; +---------------------------------------------+----------+ | ST_AsText( ST_Transform(g, 4269) ) | org_srid | +---------------------------------------------+----------+ | POINT(35.00000738369947 135.00001549175067) | 4612 | | POINT(35.00000738458632 135.00001549175067) | 4326 | | POINT(35.00325525051159 134.99723463674064) | 4301 | +---------------------------------------------+----------+
SRID:4269 というのは、NAD83 です。
mysql> SELECT SRS_NAME, SRS_ID FROM INFORMATION_SCHEMA.ST_SPATIAL_REFERENCE_SYSTEMS WHERE SRS_ID=4269; +----------+--------+ | SRS_NAME | SRS_ID | +----------+--------+ | NAD83 | 4269 | +----------+--------+
なぜ JGD2011 は変換できないのか
我々がよく使うことになりそうな JGD2011 が ST_Transform() に対応していないというのは、どういうことでしょうか。先ほどのエラーメッセージをよく見ると、「The spatial reference system has no TOWGS84 clause.」と言っています。どうやら、JGD2011には "TOWGS84" 句が記述されていないようです。
JGD2011の definition を見てみましょう*3（見やすいように整形済）
GEOGCS[
"JGD2011",
DATUM["Japanese Geodetic Datum 2011",
SPHEROID["GRS 1980",
6378137,
298.257222101,
AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","1128"]],
PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["degree",0.017453292519943278, AUTHORITY["EPSG","9122"]],
AXIS["Lat",NORTH],
AXIS["Lon",EAST],
AUTHORITY["EPSG","6668"]] 確かに、TOWGS84 なんて記述はない。。
JGD2000を見てみると、
GEOGCS[
"JGD2000",
DATUM["Japanese Geodetic Datum 2000",
SPHEROID["GRS 1980",
6378137,
298.257222101,
AUTHORITY["EPSG","7019"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6612"]],
PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["degree",0.017453292519943278, AUTHORITY["EPSG","9122"]],
AXIS["Lat",NORTH],
AXIS["Lon",EAST],
AUTHORITY["EPSG","4612"]] ある！！ある！！！ TOWGS84 が！！！
内容はオールゼロですが、これが、先ほどの WGS84 への変換実験で、緯度経度に変化がなかった理由ですね。
ちなみに、SRID:4301(Tokyo)では、
TOWGS84[-147,506,687,0,0,0,0]
となっており、TokyoからWGS84への変換では、結構変化しそうな感じが伝わってきます（7個ある数字の意味については今後調べていきたいと思います）。
なぜ JGD2011には TOWGS84 の記述がないのか
知りません。わかりません。
ST_SPATIAL_REFERENCE_SYSTEMS の値が、どこか別のところで定められた定義をそのまま持ってきているのであれば、そのおおもとがどうなっているかを辿っていく必要がありそうですし、MySQL独自に定義しているものであれば、JGD2011に TOWGS84を（オールゼロで）記述してもらうようリクエストを出せばよさそうです。
有識者からの情報に期待したいと思います :-)
TOWGS84 が記述されている座標系
mysql> SELECT SRS_NAME, SRS_ID FROM INFORMATION_SCHEMA.ST_SPATIAL_REFERENCE_SYSTEMS WHERE definition LIKE 'GEOGCS%' AND definition LIKE '%TOWGS84%' ORDER BY SRS_NAME; +-------------------------------+--------+ | SRS_NAME | SRS_ID | +-------------------------------+--------+ | Abidjan 1987 | 4143 | | Accra | 4168 | : | Yoff | 4310 | | Zanderij | 4311 | +-------------------------------+--------+ 323 rows in set (0.09 sec)
地理座標系（GEOGCS）は、MySQL 8.0.13 では、479 個が定義されているので、 323/479 に TOWGS84 が記述されているということになります。
今後の期待
JGD2011にも TO_WGS84 が記述されて、地理座標系どうしの変換仲間にいれてもらうことも勿論ですが、それ以上に 地理座標系と平面直角座標系との変換に期待したいところです。最近のMySQLのリリーススパンは結構安定してきて、その流れに従えば 次回は 2019年1月後半に MySQL 8.0.14 が出ます。ST_Transform() がどの程度進化するのか、どんなGIS関数が追加されるのか、楽しみです。機能が非常に限定された ST_Transform()関数。いま理解しておくと、MySQLの進化とともにあなたも成長します。たぶん。 いま始めれば、あなたも ST_Transform()関数のプロフェッショナル(笑)。ぜひお試しください。
※写真は「さまざまな135度」
