MySQLのJGD2011の SRS定義には towgs84が含まれていないために座標系の変換や換算できない課題について、「EPSG 9.0や9.1には、JGD2011にtowgs84の定義は含まれていない。EPSG v10.042 以降をMySQL開発チームが採用してくれれば全て解決する」との予想的結論をひとつ前のエントリで書きました。
sakaik.hateblo.jp
このエントリーでは一足先に、JGD2011でtowgs84を含んだ定義での動作を得るための方法を実験します。本稿は MySQL 8.1.0で動作確認しています。
SRSの定義の変更
MySQLのソースコードの、mysql_system_tables_data.sql ファイルにSRSの定義が納められているので、ここを書き換えてビルドするなどの方法がありそうですが、少し手順が戻りすぎの感があるので、ここでは既に稼働しているMySQL上でどうにかする方法を考えます。
MySQL には、CREATE OR REPLACE SPATIAL REFERENCE SYSTEM という素晴らしい命令があります。これを使いましょう。
現在の課題点の確認
課題(1) JGD2011からWGS84への変換がエラー
WGS84, JGD2000, JGD2011それぞれからの WGS84への変換を試みます。
データ作成
CREATE TABLE p (g geometry);
INSERT INTO p VALUES( ST_GeomFromText('POINT(35.5 135)',4612));
INSERT INTO p VALUES( ST_GeomFromText('POINT(35.5 135)',4326));
INSERT INTO p VALUES( ST_GeomFromText('POINT(35.5 135)',6668));データ確認
mysql> SELECT ST_AsText(g), ST_SRID(g) FROM p; +-----------------+------------+ | ST_AsText(g) | ST_SRID(g) | +-----------------+------------+ | POINT(35.5 135) | 4612 | | POINT(35.5 135) | 4326 | | POINT(35.5 135) | 6668 | +-----------------+------------+ 3 rows in set (0.01 sec)
これらのデータを ST_Transform()関数を使って SRID 4326 への変換を試みると、SRID 6668からの変換ができない旨エラーとなります。
mysql> SELECT ST_AsText(ST_Transform(g, 4326)), ST_SRID(ST_Transform(g, 4326)) FROM p; ERROR 3743 (22S00): Transformation from SRID 6668 is not supported. The spatial reference system has no TOWGS84 clause.
念のため、SRID 6668の行を除外して実施すれば期待通りに動作するので、SRID 6668の行が問題であることが確認できました。
mysql> SELECT ST_AsText(ST_Transform(g, 4326)), ST_SRID(g) fromsrs, ST_SRID(ST_Transform(g, 4326)) tosrs FROM p
-> WHERE ST_SRID(g)<>6668;
+----------------------------------+---------+-------+
| ST_AsText(ST_Transform(g, 4326)) | fromsrs | tosrs |
+----------------------------------+---------+-------+
| POINT(35.5 135) | 4612 | 4326 |
| POINT(35.5 135) | 4326 | 4326 |
+----------------------------------+---------+-------+
2 rows in set (0.00 sec)
課題(2) 平面直角座標系との換算がエラー
JGD2000の平面直角座標系2系(SRID 2444)から地理座標計への換算を確認します。末尾数値が少し気になりますが、正しく2系原点の緯度経度を返してくれることが確認できます。
mysql> SELECT ST_AsText(ST_Transform( ST_GeomFromText('POINT(0 0)', 2444), 4612));
+---------------------------------------------------------------------+
| ST_AsText(ST_Transform( ST_GeomFromText('POINT(0 0)', 2444), 4612)) |
+---------------------------------------------------------------------+
| POINT(33.00000000000004 131.0000000000001) |
+---------------------------------------------------------------------+
1 row in set (0.00 sec)同様にして JGD2000の平面直角座標系2系から地理座標計への換算を試みると、エラーになります。towgs84が含まれていないためです。
mysql> SELECT ST_Transform( ST_GeomFromText('POINT(0 0)', 6670), 6668);
ERROR 3743 (22S00): Transformation from SRID 6670 is not supported. The spatial reference system has no TOWGS84 clause.
JGD2011の定義を変えちゃう！
ST_SPATIAL_REFERENCE_SYSTEMSテーブルのに格納されている定義を変更してしまいます。やりたいことは、JGD2011の定義に、JGD2000にあるのと同様の形となるよう TOWGS84 の情報を追加することです。 ST_SPATIAL_REFERENCE_SYSTEMS テーブルは INFORMATION_SCHEMAにあり、一見普通のテーブルに見えますが、これを直接書き換えることはできません。専用の命令「CREATE OR REPLACE SPATIAL_REFERENCE_SYSTEM」文を使用します。
まず現在のJGD2011の定義を確認します。今回は 地理座標計と2系だけを見れば良かったのですが、ちょっと目視しておこうと、周辺のものも少し見てみました。
SELECT * FROM information_schema.ST_Spatial_reference_systems WHERE srs_id in (6668,6669,6670,6671); +------------------------------------------+--------+--------------+--------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+-------------+ | SRS_NAME | SRS_ID | ORGANIZATION | ORGANIZATION_COORDSYS_ID | DEFINITION | DESCRIPTION | +------------------------------------------+--------+--------------+--------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+-------------+ | JGD2011 | 6668 | EPSG | 6668 | 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"]] | NULL | | JGD2011 / Japan Plane Rectangular CS I | 6669 | EPSG | 6669 | PROJCS["JGD2011 / Japan Plane Rectangular CS I",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"]],PROJECTION["Transverse Mercator",AUTHORITY["EPSG","9807"]],PARAMETER["Latitude of natural origin",33,AUTHORITY["EPSG","8801"]],PARAMETER["Longitude of natural origin",129.5,AUTHORITY["EPSG","8802"]],PARAMETER["Scale factor at natural origin",0.9999,AUTHORITY["EPSG","8805"]],PARAMETER["False easting",0,AUTHORITY["EPSG","8806"]],PARAMETER["False northing",0,AUTHORITY["EPSG","8807"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",NORTH],AXIS["Y",EAST],AUTHORITY["EPSG","6669"]] | NULL | | JGD2011 / Japan Plane Rectangular CS II | 6670 | EPSG | 6670 | PROJCS["JGD2011 / Japan Plane Rectangular CS II",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"]],PROJECTION["Transverse Mercator",AUTHORITY["EPSG","9807"]],PARAMETER["Latitude of natural origin",33,AUTHORITY["EPSG","8801"]],PARAMETER["Longitude of natural origin",131,AUTHORITY["EPSG","8802"]],PARAMETER["Scale factor at natural origin",0.9999,AUTHORITY["EPSG","8805"]],PARAMETER["False easting",0,AUTHORITY["EPSG","8806"]],PARAMETER["False northing",0,AUTHORITY["EPSG","8807"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",NORTH],AXIS["Y",EAST],AUTHORITY["EPSG","6670"]] | NULL | | JGD2011 / Japan Plane Rectangular CS III | 6671 | EPSG | 6671 | PROJCS["JGD2011 / Japan Plane Rectangular CS III",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"]],PROJECTION["Transverse Mercator",AUTHORITY["EPSG","9807"]],PARAMETER["Latitude of natural origin",36,AUTHORITY["EPSG","8801"]],PARAMETER["Longitude of natural origin",132.177777777778,AUTHORITY["EPSG","8802"]],PARAMETER["Scale factor at natural origin",0.9999,AUTHORITY["EPSG","8805"]],PARAMETER["False easting",0,AUTHORITY["EPSG","8806"]],PARAMETER["False northing",0,AUTHORITY["EPSG","8807"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",NORTH],AXIS["Y",EAST],AUTHORITY["EPSG","6671"]] | NULL | +------------------------------------------+--------+--------------+--------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+-------------+
別途、JGD2000の定義を見て、TOWGS84を入れる場所を確認しておきます。
また、その測地系が使用されていると更新できないらしいので、念のため、使用されていないことも確認しておきます。
-- JGD2011(geom) mysql> SELECT * FROM INFORMATION_SCHEMA.ST_GEOMETRY_COLUMNS WHERE SRS_ID=6668; Empty set (0.00 sec) -- JGD2011(平面直角2系) mysql> SELECT * FROM INFORMATION_SCHEMA.ST_GEOMETRY_COLUMNS WHERE SRS_ID=6670; Empty set (0.00 sec)
さて、定義の書き換えです！ まず地理座標計のほう。 DEFINITIONにTOWGS84を加えた以外は元の値をそのままコピーして記述しています。WARNINGが出るのは、EPSGの予約範囲のID（0～32767と6000万台）および20億台の一部に該当するIDであるためです。warningは出るけど正常にデータは更新されます。
mysql> CREATE OR REPLACE SPATIAL REFERENCE SYSTEM 6668
-> NAME 'JGD2011'
-> ORGANIZATION 'EPSG' IDENTIFIED BY 6668
-> DEFINITION 'GEOGCS["JGD2011",DATUM["Japanese Geodetic Datum 2011",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],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"]]';
Query OK, 0 rows affected, 1 warning (0.01 sec)
mysql> show warnings;
+---------+------+--------------------------------------------------------------------------------------------------------------------------------------------------+
| Level | Code | Message |
+---------+------+--------------------------------------------------------------------------------------------------------------------------------------------------+
| Warning | 3715 | The SRID range [0, 32767] has been reserved for system use. SRSs in this range may be added, modified or removed without warning during upgrade. |
+---------+------+--------------------------------------------------------------------------------------------------------------------------------------------------+
1 row in set (0.00 sec)同様に平面直角座標2系の定義も更新。
mysql> CREATE OR REPLACE SPATIAL REFERENCE SYSTEM 6670
-> NAME 'JGD2011 / Japan Plane Rectangular CS II'
-> ORGANIZATION 'EPSG' IDENTIFIED BY 6670
-> DEFINITION 'PROJCS["JGD2011 / Japan Plane Rectangular CS II",GEOGCS["JGD2011",DATUM["Japanese Geodetic Datum 2011",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],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"]],PROJECTION["Transverse Mercator",AUTHORITY["EPSG","9807"]],PARAMETER["Latitude of natural origin",33,AUTHORITY["EPSG","8801"]],PARAMETER["Longitude of natural origin",131,AUTHORITY["EPSG","8802"]],PARAMETER["Scale factor at natural origin",0.9999,AUTHORITY["EPSG","8805"]],PARAMETER["False easting",0,AUTHORITY["EPSG","8806"]],PARAMETER["False northing",0,AUTHORITY["EPSG","8807"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",NORTH],AXIS["Y",EAST],AUTHORITY["EPSG","6670"]]';
Query OK, 0 rows affected, 1 warning (0.01 sec)
mysql> show warnings;
+---------+------+--------------------------------------------------------------------------------------------------------------------------------------------------+
| Level | Code | Message |
+---------+------+--------------------------------------------------------------------------------------------------------------------------------------------------+
| Warning | 3715 | The SRID range [0, 32767] has been reserved for system use. SRSs in this range may be added, modified or removed without warning during upgrade. |
+---------+------+--------------------------------------------------------------------------------------------------------------------------------------------------+
1 row in set (0.00 sec) これでJGD2011の地理座標と2系のSRS定義が更新されました。
いざ、JGD2011での変換／換算へ！
ST_Transform()での WGS84への変換。
mysql> SELECT ST_AsText(ST_Transform(g, 4326)), ST_SRID(g) fromsrs, ST_SRID(ST_Transform(g, 4326)) tosrs FROM p; +----------------------------------+---------+-------+ | ST_AsText(ST_Transform(g, 4326)) | fromsrs | tosrs | +----------------------------------+---------+-------+ | POINT(35.5 135) | 4612 | 4326 | | POINT(35.5 135) | 4326 | 4326 | | POINT(35.5 135) | 6668 | 4326 | +----------------------------------+---------+-------+ 3 rows in set (0.00 sec)
まったく問題ありません。
平面直角2系から地理座標への換算
JGD2011での換算もエラーなく実施できるようになりました！
mysql> SELECT ST_AsText(ST_Transform( ST_GeomFromText('POINT(0 0)', 6670), 6668));
+---------------------------------------------------------------------+
| ST_AsText(ST_Transform( ST_GeomFromText('POINT(0 0)', 6670), 6668)) |
+---------------------------------------------------------------------+
| POINT(33.00000000000004 131.0000000000001) |
+---------------------------------------------------------------------+
1 row in set (0.00 sec)
実はこの、平面直角座標系から地理座標への変換、色々おかしいぞ？という動作をしているのですが、これはもう少し調査をしてから話題にしたいと思います。ちなみに以下のようなクエリで、2つほど「おかしいぞ」があります:-)
mysql> SELECT ST_AsText(ST_Transform( ST_GeomFromText('POINT(-10000 0)', 6670), 6668));
+--------------------------------------------------------------------------+
| ST_AsText(ST_Transform( ST_GeomFromText('POINT(-10000 0)', 6670), 6668)) |
+--------------------------------------------------------------------------+
| POINT(32.99995413312006 130.8929839645656) |
+--------------------------------------------------------------------------+
まとめ
SRS定義を変更することで、JGD2011定義にも TOWGS84を加えることができることの確認を行いました。
実際の運用に於いては、この定義テーブルを書き換えることのリスク（バージョンアップ時に元に戻ってしまう可能性、それを確認して全DBサーバの定義を再度書き直す手間とその際の操作ミス等のリスクなど）が許容できる場合には、有効な手段になるのではないでしょうか。
なお、今回は平面直角座標系は2系のみを定義変更しましたが、言うまでもなく、実際に運用される場合は1～19系全ての定義に TOWGS84を加える必要があることに留意ください。
.
.
先ほどの「おかしいぞ」の答えを書いておくと
- X軸のみを南に10km移しているのに、もう一方の軸の値も変わってしまっているぞ？（実際は様々な距離で試して、値がどう変化するのかを確認しています）
- 第1引数はX軸。X軸って、、、どっちだっけ？ そう。北側ですね。変化した軸は・・・経度！？
→追記(2023/09/09)：X Y の順序はともかく、1番目の「X軸のみを南に移動」の件は、これで正しいことがわかりました。平面直角座標系に関する私の知識がまったく不足していたことが原因で、以下のエントリを読むと平面直角座標系がよくわかりますのでぜひご覧ください。
sakaik.hateblo.jp
