EPSG の towgs84 問題について

坂井さんによる以下の記事で、 MySQL に JGD2011 の定義（ WKT ）に `towgs84` が含まれない問題についてまとめられていました。

https://sakaik.hateblo.jp/entry/20230830/why_jgd2011_on_mysql_doesnd_have_towgs84_definition

MySQL に限らず、この問題は前々から一部の人に知られていましたが、改めて近年の状況を調べてみました。

自身の時系列に沿っているため、読みづらいと思います。

## 前提知識

GIS の座標系変換では、おおまかに

1. 投影座標系 A で使われている地理座標系 A' に逆投影
1. 地理座標系 A' で使われている地心直交座標系 A'' に変換
1. 地心直交座標系 A'' から WGS84 における地心直交座標系に変換（原点、軸の向きの変換）
1. WGS84 における地心直交座標系から地心直交座標系 B'' に変換（原点、軸の向きの変換）
1. 地心直交座標系 B'' から基準楕円体を設定し地理座標系 B' に変換
1. 地理座標系 B' から指定の投影法により投影座標系 B に投影

といった流れで、 WGS84 地心直交座標系の座標値を介して行われます。

人工衛星の軌道情報などによって精密な地球重心（地心）の推定値を得ることができなかった時代の測地系では、今日の WGS84 とは地心がずれていることがあります。原点がどのくらいズレており、また軸の向きはどうズレているのか。という情報を `towgs84` パラメータで示します。

`towgs84` は Helmert 変換における3パラメータもしくは7パラメータで与えられます。

https://proj.org/en/latest/operations/conversions/latlon.html#cmdoption-arg-towgs84

https://zenn.dev/boiledorange73/articles/0036-towgs84

## 今までの知識

この `towgs84` 問題は前々から存在していました。

一旦、調査する前の自身の認識を整理してみます。

### 各測地系（における地理座標系の動向）

日本でいえば、旧日本測地系 (Tokyo Datum) はベッセル楕円体が使われていたことが有名ですが、原点のズレもありました（これだけズレがあると変換後の誤差がもっとも少なくなると算出されたズレの量）。このため、 EPSG:4301 Tokyo (long/lat) では 

```:EPSG:4301 / Tokyo
+proj=longlat +ellps=bessel +towgs84=-146.414,507.337,680.507,0,0,0,0 +no_defs
```

のように `towgs84` を含んで定義されています。ズレを補正しないといけないので当然ですね。

EPSG:4612 JGD2000 (long/lat) では

```:EPSG:4612 / JGD2000
+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) には

```:EPSG:6668 / JGD2011
+proj=longlat +ellps=GRS80 +no_defs
```

のように `towgs84` のパラメータそのものが記述されていないというものでした。

### Proj の仕様

FOSS4G の座標系変換ライブラリ / ユーティリティである Proj には測地系の変換を行う際、変換元と変換先の両方に `towgs84` パラメータが含まれていない場合（たとえ `towgs84=0,0,0,0,0,0,0` に相当するものであっても）測地系変換自体を行わないとする仕様が Proj 4.6 以降ありました。

https://proj.org/en/5.0/faq.html#why-do-i-get-different-results-with-4-5-0-and-4-6-0

https://qiita.com/kochizufan/items/f4dcec7ab8e1167870ee

### そのため生じる問題

そのため、 EPSG:4301 Tokyo から直接 EPSG:6668 JGD2011 に変換（投影表示含む）しようとすると、 JGD2011 に `towgs84` のパラメータが存在しないため測地系変換が行われない。という問題がありました。

https://groups.google.com/g/qgisshitumon01/c/9afZj2IDzZY/m/mleSVPysAAAJ

## 今回の MySQL の件

それを踏まえ、今回の件を聞いたとき、同じ問題が生じているのかな？　と思って改めて調べてみると、以前とは少し状況が変化していました。

### epsg.io の表示

MapTiler Team が運営している [epsg.io](https://epsg.io/) ですが、 EPSG:4612 も EPSG:6668 も、 Proj 記法における定義では `towgs84` の記載がありませんでした。

私の記憶では EPSG:4612 は付いているけど EPSG:6668 には付いていなかったはずなのですが。。。

https://epsg.io/4612

https://epsg.io/6668

### Proj のデータベース

手元のマシンに入っている Proj 9.2.1 で確認してみると、 Proj 記法では両方とも `towgs84=0,0,0,0,0,0,0` が含まれています。

```shell-session:shell
$ proj
Rel. 9.2.1, June 1st, 2023
usage: proj [-bdeEfiIlmorsStTvVwW [args]] [+opt[=arg] ...] [file ...]
```

```shell-session:shell
$ 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]]
```

```shell-session:shell
$ 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` の有無によって測地系変換の仕様に変化はあるか確認してみます。

```shell-session:shell
$ 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 でもちゃんと変換ができるようになっています。

```shell-session:shell
$ 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](https://epsg.io/) を参考とすることが多いですが、本家である [epsg.org](https://epsg.org/) の状況を確認してみます。（いつのまにか epsg-registry.org でなくなってたんですね）

ウェブ上でパラメータなどを見ることもできますが、はっきり Proj 記法や WKT などは掲載されていません。

ユーザ登録が必要ですが、登録されているデータセットの生データを無料でダウンロードが可能です。現時点での最新版は EPSG v10.094 でしたが、 MS Access database / MySQL script / Oracle script / PostgreSQL script / WKT の各形式でダウンロードすることができます。

この中から WKT 形式をダウンロードし、どのようなデータになっているか確認してみます。

```EPSG-CRS-4301.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 上の登録状況を確認します。

```shell-session:shell
$ 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 内で座標系情報を格納しているはずだと思いますが。。。

https://github.com/OSGeo/PROJ/tree/master/data/sql

これらの情報を見ながら、ふと、座標系自体を Proj 形式や WKT 形式で定義されているのではなく、様々な定義の断片より（少なくとも Proj ではビルド時でなく）実行時に動的に生成されているのでは？　という疑念が浮かび、そっち方面に気を配りながら Proj のソースを眺めてみました。

https://github.com/OSGeo/PROJ/blob/master/src/iso19111/io.cpp

このあたりを見てると、やっぱり動的に生成している気がしてきます。

ということは Proj 記法を生成する際に、 JGD2000 には `towgs84` パラメータを生成する何らかの情報があったが、かつては JGD2011 にはそういった定義はなかったため `towgs84` が生成されなかった？

### epsg.org 巡回中

あー。これかなあ。

Transformation の定義として EPSG:1826 JGD2000 to WGS 84 (1) が登録されています。 CHANGE ID を見ると新しそうですが、 CHANGE ID 情報を上書きしてしまったのか（？）アーカイブや Proj の過去のソースを確認すると古くからちゃんと存在することが見受けられます。

https://epsg.org/transformation_1826/JGD2000-to-WGS-84-1.html

そして EPSG:9936 JGD2011 to WGS 84 (1) が2021年12月に作成されたようです。

https://epsg.org/transformation_9936/JGD2011-to-WGS-84-1.html

おそらく、 `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 は

```EPSG:4301 Tokyo (整形済み)
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)` です。

https://epsg.org/transformation_1305/Tokyo-to-WGS-84-5.html

一方で、もっともよく知られているパラメータをもつものは EPSG:15484 Tokyo to WGS 84 (108) で、 `(-146.414, 507.337, 680.507, 0.0, 0.0, 0.0, 0.0)` です。

https://epsg.org/transformation_15484/Tokyo-to-WGS-84-108.html

おそらく MySQL では (5) を用い、 Proj では (108) のバリアントを利用して `towgs84` を生成しているのだと思われます。

バリアント (108) のパラメータはもともと国土地理院の飛田幹男氏（元地理院長）の「[最近の測地座標系と座標変換についての考察](https://www.jstage.jst.go.jp/article/sokuchi1954/43/4/43_4_231/_pdf)」を出典とする EPSG:15483 Tokyo to JGD2000 (1) をもとにし、 JGD2000 と WGS84 が同等とみなしてパラメータ流用したもののようです。

#### TKY2JGD

ただし、地殻変動や地震による変動は日本列島全体で均質的に生じているわけではないため、どちらのバリアント由来の `towgs84` パラメータを用いたとしても、 Helmert 変換による測地系変換では平均的な変換でしかないため、各地点の座標を精密に変換することはできません。

高精度に変換を行うには国土地理院が提供する TKY2JGD や、それで使用されている各地点の変動量をまとめたデータ（パラメータファイル）を利用する必要があります。

https://www.gsi.go.jp/sokuchikijun/tky2jgd.html

https://www.gsi.go.jp/sokuchikijun/tky2jgd_download.html

また（データ形式は異なりますが）各地点の変動量ファイルを指定して変換する NTv2 などのグリッド変換の方法が Proj では可能です。

たとえば ESRI ArcGIS では独自で国土地理院のパラメータファイルを tky2jgd.gsb ファイルに変換し、旧日本測地系からの変換ツールを提供しているようです。その関係で、 EPSG の Transformation に tky2jgd.gsb を利用する NTv2 変換による EPSG:6712 Tokyo to JGD2000 (2) が登録されています。

https://epsg.org/transformation_6712/Tokyo-to-JGD2000-2.html

同じように、 touhokutaiheiyouoki2011.gsb を利用する EPSG:6713 JGD2000 to JGD2011 (1) が登録されています。

https://epsg.org/transformation_6713/JGD2000-to-JGD2011-1.html

ただし、 JGD2000 から JGD2011 の測地系変換では本来であれば岩手・宮城内陸地震などの JGD2000 の元期である1997.0から JGD2011 の元期の2011.4（東日本のみ）のあいだのパラメータも反映すべきところです。 ESRI が用意した touhokutaiheiyouoki2011.gsb が全部含んだものなのかは不明です。

https://www.gsi.go.jp/sokuchikijun/sokuchikijun41012.html

国土地理院のパラメータファイルを NTv2 の gsb ファイルに変換したものは下記においています。（ TKY2JGD.gsb に関しては元のパラメータファイル名に従い、大文字です。）

https://github.com/tohka/JapanGridShift

このデータに関して PROJ-data に含めるとか含めないとかいう話もあって、権利関係（測量法関係）もクリアしてる旨は示しているのですが、結局取り込まれないままよくわからない状態です。

というわけで、あくまで野良データですが、 Proj や QGIS では gsb ファイルを指定する（指定の場所に置く）ことで NTv2 による変換も可能です。

https://qiita.com/tohka383/items/e73c7235efc15efe2c1b

## まとめ

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
    - [Cape to WGS 84 (1)](https://epsg.io/1128)
- Datum
  - Geodetic
    - [Japanese Geodetic Datum 2011](https://epsg.io/1128-datum)
- Area (Extent)
  - [Jamaica](https://epsg.io/1128-area)

のように複数のカテゴリでマッチします。このため epsg.io でも epsg.org でも検索結果のページ上でカテゴリを選択することができたり、検索フィルタを設定することができたりします。

まあ、各テーブルにおける ID と考えれば重複していてもよいのですが、各文脈ごとに EPSG に登録された ID 何番。みたいな意味合いで `EPSG:xxxx` と表記するためややこしいことになっています。