【GIS】pyproj で標高の計算

はじめに

標高の計算はPROJでできるのは、下記のQiitaの記事で知っていました。

おじさんもこれは以前に確認していました。

$ echo 34.290 135.630 0.0 | PROJ_NETWORK=ON cs2cs -f "%.4f" EPSG:6667 EPSG:6697
34.2900 135.6300 -39.8601

EPSG:6667はJGD2011(Japan, Geodetic Datum 2011, 日本測地系2011)に対応するEPSGコードです。一方、EPSG:6697はJGD2011にJGD2011(Vercical)として日本の標高に対応する参照系となっています。つまり、geoid高さを楕円体高さから引き算して標高になっている、というわけです。WGS84ということならEPSG:4326を参照しても良いようです。

実はこのCUIプログラム cs2cs を直接呼ばなくても、pyproj で標高計算ができることを確認したものでメモしておきます。

PYPROJ

動作確認をとった環境

何度かこちらにも投稿していますが、座標系の変換を行ってくれるとても便利なライブラリでありプログラム群である PROJ にはpython libary として pyrpoj として公開されています。昔はPROJや pyproj のインストールに悩みましたが、最近はインストール作業をしていないせいか、昔は何が問題だったか忘れてしまいました。（多分、ライブラリのインストールとそのための環境変数の設定やパスを通すとか。。。^^;）

私のPC(Windows ノートPCでWSL2)の設定:

$ cat /etc/lsb-release
DISTRIB_ID=Ubuntu
DISTRIB_RELEASE=20.04
DISTRIB_CODENAME=focal
DISTRIB_DESCRIPTION="Ubuntu 20.04.2 LTS"
$ python -V
Python 3.8.10
$ proj
Rel. 8.1.0, July 1st, 2021
usage: proj [-bdeEfiIlmorsStTvVwW [args]] [+opt[=arg] ...] [file ...]

PROJ_NETWORK

PROJ_NETWORK_ONに対応する設定をする。現在の状態は、pyproj.network.is_network_enagled()で見れる。

import pyproj
pyproj.network.set_network_enabled(True)

環境変数を設定できるのかもしれませんが、確実なのはプログラムでset_network_enabled(True) のが良いかも。

日本のジオイドモデルで変換する

Transformer のインスタンスを作って transform するだけです。

tr = pyproj.Transformer.from_crs(6667,6669)
p = tr.transform(34.290, 135.630, 0.0)
print(p)

で以下が出力されます。

(34.29, 135.63, -39.8601222393958)

日本以外の場所を入れたらどうなるか？と気になる方もいらっしゃると思いますが、軽くスルーされるだけでした。

In []: p = tr.transform(1.0, 2.0, 3.0)

In []: p
Out[]: (1.0, 2.0, 3.0)

世界のジオイドモデルで変換する

あとは使うEPSGコードさえ分かれば、他のジオイドモデルでも変換できます。
世界をカバーするジオイドも出るEGS2008では、EPSG:6893で動作確認とれました。JGD2011(height)とは少し値が異なりますね。

In []: tr = pyproj.Transformer.from_crs(4326,6893)

In []: tr.transform(34.290, 135.630, 0.0)
Out[]: (15098262.536291692, 4043736.5764537766, -39.56243024291974)

で、それぞれのジオイドモデルとの変換ですが、とりあえず動かしてみたもののリストです。
ここではEGM2008を参照するEPSG3855のように垂直なモデルだけでなく、水平方向のモデルも含んだものを利用しています。

EPSG	地域	名前
6893	世界	EGM2008
?	イギリス	OSGB 1936 / British National Grid -- United Kingdom Ordnance Survey
6697	日本	JGD2011 + JGD2011 (vertical) height

おまけ

何とかEGM2008で変換はできたが、各国のジオイオ度を用いる計算はまだできていない。

役に立ちそうなサイト

• GEOID Calculator: UNAVOCOにあった役に立ちそうなページ。
• PROJの Vertical CRS: これを調べるべきだったのかも。理解できていない。
• EPSGの中身の調べ方：projinfo でしら得られる。例： $ projinfo epsg:6669 -o WKT2:2019
• ここ[https://qiita.com/Olivine-Ryo/items/7198ee768eb2b26492a4]によるとNETWORKについてはproj_data があれば不要になるとか。理解できていない。
• vertical datum

将来の課題として、

• 各国のGeoidoのEPSGと複合的なEPSG の整理
• GEOIDモデルを参照したいだけなのに、わざわざcompound な設定をしているCRSを参照するのは変な気がする。良い手があるのか要調査。
• EPSGの使いこなし。範囲内の判定とか。

今週も生き延びれるかな。もう疲れてきた。。。
(2021/10/10)