gmt grdmathの使い方

grdmath — GMT 6.3.0 documentation

この記事はモジュール説明とは言えないほど大雑把なものです。
深く学びたい方は公式ドキュメントをご覧ください。

大人の事情なやつです。察してください。

gmt grdmath はグリッドデータの編集・計算ができるモジュールです。
数学的なグリッドなら使えそうです。DEM等ならば、gdal_calc.pyを使った方がいい気がします。

gmt grdmath 【計算式】 = 【出力ファイル名】
計算式は逆ポーランド記法
足し算はADD、引き算はSUB、掛け算はMUL、割り算はDEV
NEGを使うと正負が反転する

逆ポーランド記法については下記を見てください

逆ポーランド記法 - Wikipedia

ここでは、DEMを使って遊んでみます。グラフ等で色々やってみたい方は下記リンクが参考になります。

GMT関係の神サイトこと「いちからはじめるGMT」

いちからはじめる GMT その64 - grdmask, grdmath

DEMを用意する

なんとなく好きなところのDEMを取得してください。
gmt grdcutを使ってもいいですが、SRTM1でも30mメッシュなので細かい方がいいでしょう。

基盤地図情報は5mメッシュなのでおすすめ。

基盤地図情報ダウンロードサービス

JPGISからGeoTIFF形式にしておくといいでしょう。
変換にはエコリスが作成している神ツールがオススメです。

基盤地図情報 標高DEMデータ変換ツール | コンテンツ | 株式会社エコリス

DEMを緯度経度にする

gmt grdcutを使った人はスルーで構いません。
基盤地図情報で作成した方はEPSG:4326に変換しておきます。

GeoTiffを加工する

GeoTiffにしたら遊んでみます。

DEMを反転する

試しに、DEMの標高値をひっくり返してみましょう。
正負を反転させるにはNEGを使います、

入手したGeoTiffを base.tif という名前にして、

gmt grdmath "base.tif" NEG = "sakasama.nc"

とすれば、標高が逆転したDEMが作成されます。

SAGAを使ってみる

SAGA GISのインストール
SAGA GISのインストール方法について紹介しています。

今回はSAGAも使ってみましょう。Morphometric Protection Indexを使ってみます。

#!/bin/bash
saga_cmd ta_morphometry 7 -DEM "base.tif"  -PROTECTION './out.tif' -RADIUS 200

とすると、よく分からない画像 out.tif ができます。

先ほどの標高を逆転させたDEMにもMorphometric Protection Indexを使ってみます。

#!/bin/bash
saga_cmd ta_morphometry 7 -DEM "sakasama.nc"  -PROTECTION './sakasama_out.tif' -RADIUS 200

とすると、こちらもよく分からない画像 sakasama_out.tif ができます。

DEMの四則演算をしてみる

out.tifsakasama_out.tifで引き算と割り算をしてみましょう。

gmt grdmath "out.tiff" "sakasama_out.tiff" SUB 0.5 MUL = nazo.nc

これは (out.tif - sakasama_out.tif) * 0.5 を意味します。
nazo.ncが出力されます。

傾斜図を作成してみる

なんとなく傾斜図を作成してみます。

gdaldem hillshade "base.tif" "hillshade.tif"

hillshade.tif が傾斜図です。

まとめ

正直、DEMの計算はgdal_calc.pyの方が良いです。
GMTで計算させると座標系が狂います。

gmt grdmath "base.tif" NEG = "sakasama.nc"

gdal_calc.py -A "base.tif" --outfile=sakasama.tif --calc="A*-1"

と書けるし、

gmt grdmath "out.tiff" "sakasama_out.tiff" SUB 0.5 MUL = nazo.nc

gdal_calc.py -A './out.tif' -B './out_sakasama.tif' --outfile=onetani.tif --calc="(A-B)/2"

と書けて、座標系も維持されるので便利です。

今回はここまで。

コメント

タイトルとURLをコピーしました