EPSG を紐解く (6)

目次

1. EPSG を紐解く (1) 〜 はじめに
2. EPSG を紐解く (2) 〜 基準系、楕円体、本初子午線
3. EPSG を紐解く (3) 〜 座標参照系（測地座標参照系、鉛直座標参照系）、座標系
4. EPSG を紐解く (4) 〜 座標参照系（投影座標参照系）、投影演算
5. EPSG を紐解く (5) 〜 座標参照系（複合座標参照系）
6. EPSG を紐解く (6) 〜 座標変換、 towgs84 句

基準系間の変換 (Transformation)

たとえば日本測地系2000と日本測地系2011間の変換や、北米の NAD27 から NAD83 への変換など異なる基準系に変換を行うことがあります。 EPSG ではそのような演算を変換 (Transformation) と呼び、換算 (Conversion) と区別しています。

換算 (Conversion) では地理座標参照系などのベースとなる座標参照系があり、それをもとに換算することで投影座標参照系などの派生の座標参照系が作成されると考えることもでき、その作成された派生座標参照系は換算情報とともに個別の座標参照系として登録することができます。

一方で、変換 (Transformation) は基準系間の関係性を表しており、特定の座標参照系に付随するものではありません。また、変換方法は1種類だけとは限らず、異なる変換方法や別の変換パラメタによるバリアントが登録されていることもあります。

たとえば下記は CRS:4612 JGD2000 (geog2D) から CRS:6668 JGD2011 (geog2D) への変換になります。

source_crs_code : 4612, target_crs_code : 6668 をもつ変換は Operation:6713 JGD2000 to JGD2011 (1) と Operation:6698 JGD2000 to JGD2011 (2) の2つのバリアントが登録されています。

Operation:6713 の方は Method:9615 NTv2 すなわち別途基準系間の差異を格納したグリッドファイル (touhokutaiheiyouoki2011.gsb) を用いる方法であり、局所的で不規則な差異に対し変換が可能で、精度が 0.2m となっています。 ESRI 社の ArcGIS の変換ツールで使われている関係で登録されたとみられます（実際には東日本大震災だけでなく他の地震も加味する必要があり、不十分な変換）。

一方 Operation:6698 の方は Method:9603 Geocentric translations (geog2D domain) すなわち JGD2000 の地理座標参照系 →　JGD2000 の地心座標参照系 →　地心座標参照系の原点のシフト（シフト量はゼロ） →　JGD2011 の地心座標参照系 → JGD2011 の地理座標参照系、という順で変換する方法で精度は 1.0m とされています。

それぞれの変換で用いられるシフト量といったパラメタは換算のときと同様に epsg_coordoperationparamvalue テーブルに格納されています。

towgs84 句

変換は各座標参照系に付随するものではないと前述しました。しかし変換処理を行うたびに複数存在する変換方法を解決するのではなく、変換方法の情報を座標参照系にあらかじめ紐付けた形式もあり、そのような形式の座標参照系を Bound CRS と呼びます。もっとも代表的なものに、 towgs84 があります。これは WGS84 座標参照系をハブとなる座標参照系とし、 WGS84 座標参照系への変換方法を3パラメタまたは7パラメタで表したものです。 towgs84 を用いることで、互いの基準系間の変換が登録されていなくても変換することが可能となります。

たとえば JGD2000 の場合、 JGD2000 to JGD2011 と同様に JGD2000 to WGS84 への変換が Operation:1826 JGD2000 to WGS 84 (1) として登録されています。この変換の各シフト量は JGD2000 to JGD2011 (2) に準じて tX = 0, tY = 0, tZ = 0 であるため towgs84=0,0,0,0,0,0,0 となります。

なお JGD2011 から WGS84 への変換が最近まで登録されていなかったため、 JGD2011 に対し towgs84 句を付与した Bound CRS を作成することができませんでした。加えて PROJ では towgs84 句などの WGS84 への変換手法が明示されていないと異なる基準系間の変換が行われないという仕様があるため、旧日本測地系から日本測地系2011への変換が行われないという問題がありました。

ちなみに MySQL では使用している EPSG dataset が古いため、現時点でも towgs84 句が存在しておらず課題となっています。

バインドするバリアントの選択

JGD2000 から WGS84 の場合、 Geocentric translations の手法のバリアントは1種類でしたが、同じ変換手法のバリアントが多数あったらどうでしょうか。旧日本測地系 (Tokyo Datum) の場合がまさにそれです。

CRS:4301 Tokyo から WGS84 への変換は6つ登録されています（廃止フラグが立っているものを除く）。

6つのバリアントに対し epsg_coordoperation / epsg_coordoperationparamvalue / epsg_usage / epsg_extent 等のテーブルの連結させ重要な情報のみを抜粋したものが下表となります。¹

もっともよく使われているものはバリアント108の Operation:15484 Tokyo to WGS 84 (108) になります。この towgs84 のパラメタは Operation:15483 Tokyo to JGD2000 (1) からの流用となり、大元は飛田幹男.　最近の測地座標系と座標変換についての考察からとなるようです。

実際 PROJ では +towgs84=-146.414,507.337,680.507,0,0,0,0 が選択されます。

$ projinfo --boundcrs-to-wgs84 -o PROJ EPSG:4301
PROJ.4 string:
+proj=longlat +ellps=bessel +towgs84=-146.414,507.337,680.507,0,0,0,0 +no_defs +type=crs

$ projinfo --boundcrs-to-wgs84 -o WKT1:GDAL EPSG:4301
WKT1:GDAL string:
GEOGCS["Tokyo",
    DATUM["Tokyo",
        SPHEROID["Bessel 1841",6377397.155,299.1528128,
            AUTHORITY["EPSG","7004"]],
        TOWGS84[-146.414,507.337,680.507,0,0,0,0],
        AUTHORITY["EPSG","6301"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4301"]]

先ほどの表を見ていただければ気付くと思いますが、バリアント5の精度 (accuracy) の方がよいバリアントです。実際、 MySQL では精度の値を優先しているとみられ、バリアント5をバインドした towgs84 パラエタを採用しています。

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"]]

しかし更によく見ると、バリアント5である Operation:1305 Tokyo to WGS 84 (5) の適用範囲は Extent:3266 Korea, Republic of (South Korea) - onshore すなわち韓国に限定されています。歴史的経緯はあまりわかりませんが、旧日本測地系は韓国でも使われていたそうで、このバリアントは韓国に特化したものようです。

一方で PROJ では適用範囲の面積も加味してバインドする変換を決定しているようです。

標高体系の変換

JGD2024 が未対応の現時点において、 EPSG に登録されている JGD2011 関係の標高に関する変換は Operation:6699 JGD2000 (vertical) height to JGD2011 (vertical) height (1) のみです。これは JGD2000 の標高値から JGD2011 の標高値への変換です。ただし、この変換のオフセット量はゼロなので、実質無変換となります。

EPSG dataset に登録されているものは上記のみですが、 PROJ の場合は独自に PROJ:EPSG_6667_TO_EPSG_6695 という識別子で CRS:6667 JGD2011 (geog3d) から CRS:6695 JGD2011 (vertical) height の変換、すなわち JGD2011 の緯度、経度、楕円体高から標高への変換が登録されています。

この変換にはジオイドデータである jp_gsi_gsigeo2011.tif が使用されます。

$ projinfo -s EPSG:6667 -t EPSG:6695
Candidate operations found: 1
-------------------------------------
Operation No. 1:

PROJ:EPSG_6667_TO_EPSG_6695, JDG2011 to JGD2011 height, unknown accuracy, Japan - onshore mainland - Hokkaido, Honshu, Shikoku, Kyushu.

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +inv +proj=vgridshift +grids=jp_gsi_gsigeo2011.tif +multiplier=1
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

WKT2:2019 string:
COORDINATEOPERATION["JDG2011 to JGD2011 height",
    SOURCECRS[
        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,3],
                AXIS["geodetic latitude (Lat)",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["geodetic longitude (Lon)",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["ellipsoidal height (h)",up,
                    ORDER[3],
                    LENGTHUNIT["metre",1]],
            ID["EPSG",6667]]],
    TARGETCRS[
        VERTCRS["JGD2011 (vertical) height",
            VDATUM["Japanese Geodetic Datum 2011 (vertical)"],
            CS[vertical,1],
                AXIS["gravity-related height (H)",up,
                    LENGTHUNIT["metre",1]],
            ID["EPSG",6695]]],
    METHOD["Geographic3D to GravityRelatedHeight (gtx)",
        ID["EPSG",9665]],
    PARAMETERFILE["Geoid (height correction) model file","jp_gsi_gsigeo2011.tif"],
    USAGE[
        SCOPE["Not known."],
        AREA["Japan - onshore mainland - Hokkaido, Honshu, Shikoku, Kyushu."],
        BBOX[30.94,129.3,45.54,145.87]],
    ID["PROJ","EPSG_6667_TO_EPSG_6695"]]

EPSG への JGD2024 追加はまだ行われていませんが、つい先日（4月10日）、 PROJ の方では「ジオイド2024日本とその周辺」を用いた変換がコミットされました。

ここでは PROJ:EPSG_6667_TO_EPSG_6695_2024 という識別子で、ジオイドモデル jp_gsi_jpgeo2024.tif を用いた変換です。座標参照系上は JGD2024 ではなく、あくまで JGD2011 (geog3D) から JGD2011 (vertical) への変換として定義されています。

PROJ で使用可能な形式のジオイドデータは下記から入手することができます。

最後に

2025年4月中旬現在、 EPSG への JGD2024 登録はまだ行われていません。

水平座標も含め、今後どのような形で JGD2024 の対応がなされるかわかりませんが、現在の EPSG の登録体系を知ることで多少の想定が可能となるのではないでしょうか。

併せて個人的には前々からじっくり調べたいと思っていた EPSG の構造をちゃんと学ぶことができてよかったです。

1. area は $(\lambda_2-\lambda_1) (\sin\varphi_2 - \sin\varphi_1)$ で計算される面積の近似値。面積の大小を比較ができればよいため、平方メートル単位にするのに必要な地球半径の2乗を乗じたり、経度をラジアンになおしていない。 ↩