坂井さんによる以下の記事で、 MySQL に JGD2011 の定義( WKT )に towgs84
が含まれない問題についてまとめられていました。
MySQL に限らず、この問題は前々から一部の人に知られていましたが、改めて近年の状況を調べてみました。
自身の時系列に沿っているため、読みづらいと思います。
前提知識
GIS の座標系変換では、おおまかに
- 投影座標系 A で使われている地理座標系 A' に逆投影
- 地理座標系 A' で使われている地心直交座標系 A'' に変換
- 地心直交座標系 A'' から WGS84 における地心直交座標系に変換(原点、軸の向きの変換)
- WGS84 における地心直交座標系から地心直交座標系 B'' に変換(原点、軸の向きの変換)
- 地心直交座標系 B'' から基準楕円体を設定し地理座標系 B' に変換
- 地理座標系 B' から指定の投影法により投影座標系 B に投影
といった流れで、 WGS84 地心直交座標系の座標値を介して行われます。
人工衛星の軌道情報などによって精密な地球重心(地心)の推定値を得ることができなかった時代の測地系では、今日の WGS84 とは地心がずれていることがあります。原点がどのくらいズレており、また軸の向きはどうズレているのか。という情報を towgs84
パラメータで示します。
towgs84
は Helmert 変換における3パラメータもしくは7パラメータで与えられます。
今までの知識
この towgs84
問題は前々から存在していました。
一旦、調査する前の自身の認識を整理してみます。
各測地系(における地理座標系の動向)
日本でいえば、旧日本測地系 (Tokyo Datum) はベッセル楕円体が使われていたことが有名ですが、原点のズレもありました(これだけズレがあると変換後の誤差がもっとも少なくなると算出されたズレの量)。このため、 EPSG:4301 Tokyo (long/lat) では
+proj=longlat +ellps=bessel +towgs84=-146.414,507.337,680.507,0,0,0,0 +no_defs
のように towgs84
を含んで定義されています。ズレを補正しないといけないので当然ですね。
EPSG:4612 JGD2000 (long/lat) では
+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
のように towgs84=0,0,0,0,0,0,0
という記述が含まれています。
そしてなぜか EPSG:6668 JGD2011 (long/lat) には
+proj=longlat +ellps=GRS80 +no_defs
のように towgs84
のパラメータそのものが記述されていないというものでした。
Proj の仕様
FOSS4G の座標系変換ライブラリ / ユーティリティである Proj には測地系の変換を行う際、変換元と変換先の両方に towgs84
パラメータが含まれていない場合(たとえ towgs84=0,0,0,0,0,0,0
に相当するものであっても)測地系変換自体を行わないとする仕様が Proj 4.6 以降ありました。
そのため生じる問題
そのため、 EPSG:4301 Tokyo から直接 EPSG:6668 JGD2011 に変換(投影表示含む)しようとすると、 JGD2011 に towgs84
のパラメータが存在しないため測地系変換が行われない。という問題がありました。
今回の MySQL の件
それを踏まえ、今回の件を聞いたとき、同じ問題が生じているのかな? と思って改めて調べてみると、以前とは少し状況が変化していました。
epsg.io の表示
MapTiler Team が運営している epsg.io ですが、 EPSG:4612 も EPSG:6668 も、 Proj 記法における定義では towgs84
の記載がありませんでした。
私の記憶では EPSG:4612 は付いているけど EPSG:6668 には付いていなかったはずなのですが。。。
Proj のデータベース
手元のマシンに入っている Proj 9.2.1 で確認してみると、 Proj 記法では両方とも towgs84=0,0,0,0,0,0,0
が含まれています。
$ proj
Rel. 9.2.1, June 1st, 2023
usage: proj [-bdeEfiIlmorsStTvVwW [args]] [+opt[=arg] ...] [file ...]
$ projinfo EPSG:4612
PROJ.4 string:
+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs +type=crs
WKT2:2019 string:
GEOGCRS["JGD2000",
DATUM["Japanese Geodetic Datum 2000",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["Japan - onshore and offshore."],
BBOX[17.09,122.38,46.05,157.65]],
ID["EPSG",4612]]
$ projinfo EPSG:6668
PROJ.4 string:
+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs +type=crs
WKT2:2019 string:
GEOGCRS["JGD2011",
DATUM["Japanese Geodetic Datum 2011",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["Japan - onshore and offshore."],
BBOX[17.09,122.38,46.05,157.65]],
ID["EPSG",6668]]
って、え? EPSG:6668 にも含まれているの???
WKT の定義の方には両方とも記述はありませんでした。
Proj における動作の確認
Proj の座標系変換ツール cs2cs
を用いて、 towgs84
の有無によって測地系変換の仕様に変化はあるか確認してみます。
$ echo 135 35 | cs2cs +proj=longlat +ellps=bessel +towgs84=-146.414,507.337,680.507,0,0,0,0 +no_defs +to +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
134d59'49.935"E 35d0'11.51"N 0.000
$ echo 135 35 | cs2cs +proj=longlat +ellps=bessel +towgs84=-146.414,507.337,680.507,0,0,0,0 +no_defs +to +proj=longlat +ellps=GRS80 +no_defs
135dE 35dN 0.000
やはり towgs84
がないと、地心のズレはもちろん楕円体間の変換も行ってくれない仕様は変わっていないようです。
ちなみに EPSG コードを直接指定すると、 EPSG:4612 だけでなく EPSG:6668 でもちゃんと変換ができるようになっています。
$ echo 135 35 | cs2cs +init=EPSG:4301 +to +init=EPSG:4612
134d59'49.935"E 35d0'11.51"N 0.000
$ echo 135 35 | cs2cs +init=EPSG:4301 +to +init=EPSG:6668
134d59'49.935"E 35d0'11.51"N 0.000
epsg.org の登録状況
普段は軽くて便利な epsg.io を参考とすることが多いですが、本家である epsg.org の状況を確認してみます。(いつのまにか epsg-registry.org でなくなってたんですね)
ウェブ上でパラメータなどを見ることもできますが、はっきり Proj 記法や WKT などは掲載されていません。
ユーザ登録が必要ですが、登録されているデータセットの生データを無料でダウンロードが可能です。現時点での最新版は EPSG v10.094 でしたが、 MS Access database / MySQL script / Oracle script / PostgreSQL script / WKT の各形式でダウンロードすることができます。
この中から WKT 形式をダウンロードし、どのようなデータになっているか確認してみます。
GEOGCRS["Tokyo",
DATUM["Tokyo",
ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
LENGTHUNIT["metre",1,ID["EPSG",9001]],
ID["EPSG",7004]],
ID["EPSG",6301]],
CS[ellipsoidal,2,ID["EPSG",6422]],
AXIS["Geodetic latitude (Lat)",north],
AXIS["Geodetic longitude (Lon)",east],
ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]],
ID["EPSG",4301]]
EPSG:4301 の WKT の記述を見てみると towgs84
の記載がありませんでした。
ちなみにと、 EPSG:4301 の Proj 上の登録状況を確認します。
$ projinfo EPSG:4301
PROJ.4 string:
+proj=longlat +ellps=bessel +towgs84=-146.414,507.337,680.507,0,0,0,0 +no_defs +type=crs
WKT2:2019 string:
GEOGCRS["Tokyo",
DATUM["Tokyo",
ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Geodesy."],
AREA["Japan - onshore; North Korea - onshore; South Korea - onshore."],
BBOX[20.37,122.83,45.54,154.05]],
ID["EPSG",4301]]
Proj 記法では towgs84
はありますが、 WKT にはないですね。
配布されている EPSG のデータセットのうち、他の SQL script 形式のデータをダウンロードしてきてみても、 Proj 記法や WKT っぽい記述は含まれておらず、どこで定義されているのかわかりませんでした。
proj.db を確認
Proj に登録されている座標系等の情報は SQLite3 形式で proj.db に保存されています。
また proj.db を作成するためのソースは以下のような sql ファイルとなっています。 EPSG 生データセットのことはおいとおいたとしても、 projinfo
で WKT などが表示されるからには Proj 内で座標系情報を格納しているはずだと思いますが。。。
これらの情報を見ながら、ふと、座標系自体を Proj 形式や WKT 形式で定義されているのではなく、様々な定義の断片より(少なくとも Proj ではビルド時でなく)実行時に動的に生成されているのでは? という疑念が浮かび、そっち方面に気を配りながら Proj のソースを眺めてみました。
このあたりを見てると、やっぱり動的に生成している気がしてきます。
ということは Proj 記法を生成する際に、 JGD2000 には towgs84
パラメータを生成する何らかの情報があったが、かつては JGD2011 にはそういった定義はなかったため towgs84
が生成されなかった?
epsg.org 巡回中
あー。これかなあ。
Transformation の定義として EPSG:1826 JGD2000 to WGS 84 (1) が登録されています。 CHANGE ID を見ると新しそうですが、 CHANGE ID 情報を上書きしてしまったのか(?)アーカイブや Proj の過去のソースを確認すると古くからちゃんと存在することが見受けられます。
そして EPSG:9936 JGD2011 to WGS 84 (1) が2021年12月に作成されたようです。
おそらく、 towgs84
は測地系 (datum) で定義されているわけでもなく、地理座標系ででもなく、 WGS84 への変換方法 (operation / transformation) の登録情報から付加されたものだと推測されます。
そして、この JGD2011 to WGS84 の Transformation が登録される前は EPSG:6668 JGD2011 に towgs84
が存在せず、登録された現在では towgs84=0,0,0,0,0,0,0
が付加されているのだと思われます。
Tokyo to WGS 84 の確認
一方、 Tokyo の towgs84
についてみてみます。
MySQL の ST_SPATIAL_REFERENCE_SYSTEMS
に登録されている EPSG:4301 Tokyo の WKT は
GEOGCS["Tokyo",
DATUM["Tokyo",
SPHEROID["Bessel 1841",6377397.155,299.1528128,
AUTHORITY["EPSG","7004"]],
TOWGS84[-147,506,687,0,0,0,0],
AUTHORITY["EPSG","6301"]],
PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["degree",0.017453292519943278,AUTHORITY["EPSG","9122"]],
AXIS["Lat",NORTH],
AXIS["Lon",EAST],
AUTHORITY["EPSG","4301"]]
と、 TOWGS84[-147,506,687,0,0,0,0]
が含まれているようです。
Helmert 変換を行う Tokyo to WGS 84 の変換には6つの VARIANT が登録されています。それぞれには、(ちゃんと評価した結果なのかは不明ですが)精度情報も登録されており、もっとも精度がよいとされているものは EPSG:1305 Tokyo to WGS 84 (5) になっています。このときのパラメータは (-147.0, 506.0, 687.0, 0.0, 0.0, 0.0, 0.0)
です。
一方で、もっともよく知られているパラメータをもつものは EPSG:15484 Tokyo to WGS 84 (108) で、 (-146.414, 507.337, 680.507, 0.0, 0.0, 0.0, 0.0)
です。
おそらく MySQL では (5) を用い、 Proj では (108) のバリアントを利用して towgs84
を生成しているのだと思われます。
バリアント (108) のパラメータはもともと国土地理院の飛田幹男氏(元地理院長)の「最近の測地座標系と座標変換についての考察」を出典とする EPSG:15483 Tokyo to JGD2000 (1) をもとにし、 JGD2000 と WGS84 が同等とみなしてパラメータ流用したもののようです。
TKY2JGD
ただし、地殻変動や地震による変動は日本列島全体で均質的に生じているわけではないため、どちらのバリアント由来の towgs84
パラメータを用いたとしても、 Helmert 変換による測地系変換では平均的な変換でしかないため、各地点の座標を精密に変換することはできません。
高精度に変換を行うには国土地理院が提供する TKY2JGD や、それで使用されている各地点の変動量をまとめたデータ(パラメータファイル)を利用する必要があります。
また(データ形式は異なりますが)各地点の変動量ファイルを指定して変換する NTv2 などのグリッド変換の方法が Proj では可能です。
たとえば ESRI ArcGIS では独自で国土地理院のパラメータファイルを tky2jgd.gsb ファイルに変換し、旧日本測地系からの変換ツールを提供しているようです。その関係で、 EPSG の Transformation に tky2jgd.gsb を利用する NTv2 変換による EPSG:6712 Tokyo to JGD2000 (2) が登録されています。
同じように、 touhokutaiheiyouoki2011.gsb を利用する EPSG:6713 JGD2000 to JGD2011 (1) が登録されています。
ただし、 JGD2000 から JGD2011 の測地系変換では本来であれば岩手・宮城内陸地震などの JGD2000 の元期である1997.0から JGD2011 の元期の2011.4(東日本のみ)のあいだのパラメータも反映すべきところです。 ESRI が用意した touhokutaiheiyouoki2011.gsb が全部含んだものなのかは不明です。
国土地理院のパラメータファイルを NTv2 の gsb ファイルに変換したものは下記においています。( TKY2JGD.gsb に関しては元のパラメータファイル名に従い、大文字です。)
このデータに関して PROJ-data に含めるとか含めないとかいう話もあって、権利関係(測量法関係)もクリアしてる旨は示しているのですが、結局取り込まれないままよくわからない状態です。
というわけで、あくまで野良データですが、 Proj や QGIS では gsb ファイルを指定する(指定の場所に置く)ことで NTv2 による変換も可能です。
まとめ
EPSG (本家の)データセットがどうなっているのかちゃんと見たことなかったので、いろいろ知らないことばかりでした。
まとめとして
- どうやら EPSG 本家では地理座標系や投影座標系そのものの具体的な定義を Proj 記法や WKT でしているわけではない模様
- 軸 (axis) なり楕円体 (ellipsoid) なり測地系 (geodetic_datum) なり測地系変換 (helmert_transformation) なりの部品部品がちゃんと定義されており、それを組み合わせてビルド時あるいは実行時に生成しているよう
- RDB ベースの定義がなされており、 Proj でも proj.db の状態のまま利用している
- ソフトごとに記載方法が異なるのはそのためと思われる
- いつの間にか JGD2011 に
towgs84
がついてた- おそらく、それは JGD2011 to WGS 84 が定義されたため
- epsg.io で
+towgs84=0,0,0,0,0,0,0
が表示されないのは、 Proj 記法の生成処理に( Proj の変換処理の仕様を正するのなら)不備があるから? - MySQL の JGD2011 定義に不備があるのは、 WKT 生成エンジンの問題か、取り込んだ EPSG データセットが古いから?
おまけ
このように EPSG といえば地理座標系もしくは投影座標系というイメージがありますが、実際には様々な定義がなされており、それぞれでコード (ID) が付いています。更にややこしいことに数字は分けられておらず重複していることもあります。
たとえば、 EPSG:1128 であれば
- Operation
- Transformation
- Datum
- Geodetic
- Area (Extent)
のように複数のカテゴリでマッチします。このため epsg.io でも epsg.org でも検索結果のページ上でカテゴリを選択することができたり、検索フィルタを設定することができたりします。
まあ、各テーブルにおける ID と考えれば重複していてもよいのですが、各文脈ごとに EPSG に登録された ID 何番。みたいな意味合いで EPSG:xxxx
と表記するためややこしいことになっています。