datumgrid を使った座標変換

9603 ワード
QGIS foss4g Proj.4 GDAL GDAL テキストリンク

FOSS4G Advent Calendar 2019 12日目です。

座標変換するときに。。。

日本測地系2000 (JGD2000) と日本測地系2011 (JGD2011) は、基準とする回転楕円体はいずれも GRS80 であり、 GIS の定義上は同じ座標系です。このため、 GIS 上で投影変換（レイヤ、プロジェクト CRS の設定あるいは新規保存による座標系変換）を行っても座標値として変化はありません。

しかし実際には JGD2000 と JGD2011 は元期（基準となる日時）が異なります1。長期的な地殻変動はもちろん大地震によって日本列島自体が動くため、同じ地点（地物）でも全地球的に見ると位置は変化しています。つまり地図の座標軸などは変化していなくても、年間1cmの速度で動くたかし君は2000年時点と2011年時点で違う場所います。

この元期による差は（大地震の影響などもあり）日本全国で不均一であるため、全国同一の計算式で変換することはできません。このため各地点ごとの補正量がまとめられたパラメータファイルを使って、対象地点の補正を行う必要があります。
そのためのパラメータファイルや変換ツールが国土地理院から提供されています。

さて。 FOSS4G で座標系変換といえば Proj 先輩ですが、実はパラメータファイルを利用した変換も可能で、各国のパラメータファイルのリポジトリも存在しています。ただ、現時点で日本のパラメータファイルは含まれていません。（ Proj を使用している E 社の A には同梱されています。）

今回 Proj で使用できるパラメータファイルを作成し、座標の変換をしてみました。

*.par の仕様について

国土地理院のパラメータファイルは独自形式の *.par ファイルで提供されており、下記のようなテキスト形式です。

JGD2000-TokyoDatum Ver.2.1.2
MeshCode   dB(sec)   dL(sec)
46303582  12.79799  -8.13354
46303583  12.79879  -8.13749
　　：
　　：

地域メッシュのコードと緯度経度の補正値 dB, dL が秒単位で記載されています。このパラメータファイル内では、メッシュコードは長方形のポリゴンそのものを意味しているのではなく、グリッドの南西角の点の座標値を意味します。

つまり下記の行は、北緯30.9833度　東経130.6500度の地点における緯度の補正量が12.79799秒、経度の補正量が-8.13354秒であることを意味しています。

46303582  12.79799  -8.13354

NTv2 (*.gsb) の仕様について

Proj で使用可能な（水平方向の）パラメータファイルは NTv1 (*.dat), NTv2 (*.gsb), CTable2 (*.ct2) が利用できますが、主に NTv2 形式が使用されています。（正確には NTv2 という補正方法で使用するためのファイル形式。）

NTv2 はもともと北米の座標系である NAD27 から NAD83 への変換を行うためカナダで開発された変換方法及びパラメータファイルのようです。

仕様は下記のドキュメントで確認することができますが、カナダで開発されたためか経度に関しては西がプラス、東がマイナスです。また記載する補正量の順は南東から南西に向かい、最終的に北西端となる順です。以上二点を押さえればそんなに難しくありません。

*.gsb ファイルは規則正しく並んだ点の変化量を保持した形式なので、いわばラスタファイルとみることができ、 QGIS で開くこともできますし GDAL にも対応しています。

また ESRI 社からオープンソースでテキスト形式から *.gsb に変換できるツールが公開されています。

とはいえ *.gsb ファイルの仕様は非常に単純ですので、どうせ国土地理院のパラメータファイルから変換するスクリプトを書くのなら、直接 *.gsb ファイルを作成してしまった方が楽です。

変換

国土地理院の *.par ファイルを *.gsb に変換するスクリプトを書いてみました。（ TKY2JGD.par と touhokutaiheiyouoki2011.par が必要）

Proj での使い方

$ echo 135 35 | cs2cs -f "%.9f" +proj=latlong +ellps=bessel +towgs84=-146.414,507.337,680.507 +nadgrids=tky2jgd.gsb,null +to +proj=latlong +ellps=GRS80 +towgs84=0,0,0

134.997229989   35.003215692 0.000000000

+nadgrids オプションで水平グリッドシフトファイルを指定します。 ,null をつけることにより指定したグリッドファイルの範囲外でもエラーを抑制することができます。

ただ、事前に *.gsb ファイルを適切な場所に入れておかないと、当然 Proj パイセンはファイルをみつけることができないのでエラーとなります。どこに入れればいいかは、以下のコマンドのようにデバッグレベルをあげることで確認できます。

$ echo 135 35 | PROJ_DEBUG=3 cs2cs -f "%.9f" +proj=latlong +ellps=bessel +towgs84=-146.414,507.337,680.507 +nadgrids=tky2jgd.gsb,null +to +proj=latlong +ellps=GRS80 +towgs84=0,0,0 2>&1 | fgrep gsb

pj_open_lib(tky2jgd.gsb,null): call fopen(/usr/local/Cellar/proj/6.0.0/share/proj/tky2jgd.gsb,null) - failed
    grids=tky2jgd.gsb,null
pj_open_lib(tky2jgd.gsb): call fopen(/usr/local/Cellar/proj/6.0.0/share/proj/tky2jgd.gsb) - failed

今回確認した環境では /usr/local/Cellar/proj/6.0.0/share/proj に入れるとよいことがわかります。ドキュメントを軽く読む感じだと PROJ_LIB 環境変数を設定すればいいような気もしますが、うまくいきませんでした。

もしくは +nadgrids=./tky2jgd.gsb のようにファイル名ではなくパスを指定することもできます。

GDAL/OGR での使い方

たとえば、このように使います。

$ ogr2ogr -f 'ESRI Shapefile' -s_srs '+proj=latlong +ellps=bessel +towgs84=-146.414,507.337,680.507 +nadgrids=./tky2jgd.gsb,null' -t_srs '+proj=latlong +ellps=GRS80 +towgs84=0,0,0' JGD2000.shp Tokyo.shp

-s_srs で入力元の座標系を上書き（仮に *.prj ファイル等の情報があっても無視）し、 -t_srs で出力先の座標系に変換しています。

GDAL/OGR の場合もパラメータファイルを入れる場所がわかりませんでした。先程の Proj に場所に入れてもだめでした。私がテストした環境は Mac 用の KyngChaos さんのパッケージで GDAL とか入れていつつ、更に homebrew でもいれたり、よくわかんないことになっているせいかもしれません。

とりあえず +nadgrids=./tky2jgd.gsb,null のようなパス表記だと大丈夫でした。

QGIS での使い方

事前準備として、メニューの「設定」→「カスタム投影法」で新規投影法を定義します。名前はなんでもよいのですが、たとえば Tokyo (gridshift) とし、パラメータは +proj=longlat +ellps=bessel +towgs84=-146.414,507.337,680.507 +nadgrids=/path/to/tky2jgd.gsb,null とします。

例によってどこにパラメータファイルを入れるかわかりませんでしたが、とりあえずフルパス書いたら大丈夫でした。 Windows 版の OSGeo4W だとすんなりいけたかと思います。

旧日本測地系 (Tokyo Datum) のデータ（レイヤ）に対し、レイヤプロパティを開き、「ソース」タブ→「座標参照系の設定」の欄に、定義した座標系を設定します。これにより tky2jgd.gsb を利用した Tokyo (gridshift) におけるデータとして扱われます。

プロジェクトの CRS を適宜目的の CRS に変えるなり、レイヤ右クリックメニューの「エクスポート」→「地物の保存」で JGD2000 などの座標系を指定して保存します。

注意事項

ドキュメントを見る限り NTv2 の補間方法はバイリニア補間によるもので、 Proj や国土地理院の変換ツール（ TKY2JGD / PatchJGD ）とおそらく同様の手法で計算しているものだと思われますが、座標値が正確に（国土地理院の変換ツールと同精度で）変換できているか保証するものではありません。

特に公共測量の成果に対する変換を行う際には公共測量成果改定マニュアルに従い、正しく補正されているか点検を行い、計画機関の承認を得る必要があります。

今後について

国土地理院に確認したところ、パラメータファイルは基本測量の成果ではない（として取り扱わない）ので測量法における複製申請は不要であり、 CC-BY 4.0 と互換性のある国土地理院コンテンツ利用規約を遵守すればよい。と回答をいただきました。

精度まわりなど、もう少しテストを行って問題ないようならば、 Proj の datumgrid のリポジトリにコミットしようと考えています。リポジトリに正式に含まれれば、おそらく QGIS などとともに同梱されることになるかと思います。

また高さに関しては GTX という形式でファイルを作成する必要があるのですが、まだまったく調べていません。 JGD2000 から JGD2011 の補正量ももちろんですが、ジオイド高のパラメータを GTX にすることで標高の取扱いが楽に行うことができ、恩恵は大きいのではないかと考えています。（ジオイド高はおそらく測量成果扱い）

  1. 東北から岐阜長野あたりまでは東日本大震災の影響で元期が2011年5月となっていますが、それ以外の地域は JGD2000 と同じ元期である1997年1月です。 ↩