この記事はモジュール説明とは言えないほど大雑把なものです。
深く学びたい方は公式ドキュメントをご覧ください。
大人の事情なやつです。察してください。
gmt grdmath はグリッドデータの編集・計算ができるモジュールです。
数学的なグリッドなら使えそうです。DEM等ならば、gdal_calc.pyを使った方がいい気がします。
gmt grdmath 【計算式】 = 【出力ファイル名】
計算式は逆ポーランド記法。
足し算はADD、引き算はSUB、掛け算はMUL、割り算はDEV
NEGを使うと正負が反転する
逆ポーランド記法については下記を見てください
ここでは、DEMを使って遊んでみます。グラフ等で色々やってみたい方は下記リンクが参考になります。
GMT関係の神サイトこと「いちからはじめるGMT」
DEMを用意する
なんとなく好きなところのDEMを取得してください。
gmt grdcutを使ってもいいですが、SRTM1でも30mメッシュなので細かい方がいいでしょう。
基盤地図情報は5mメッシュなのでおすすめ。
JPGISからGeoTIFF形式にしておくといいでしょう。
変換にはエコリスが作成している神ツールがオススメです。
DEMを緯度経度にする
gmt grdcutを使った人はスルーで構いません。
基盤地図情報で作成した方はEPSG:4326に変換しておきます。
GeoTiffを加工する
GeoTiffにしたら遊んでみます。
DEMを反転する
試しに、DEMの標高値をひっくり返してみましょう。
正負を反転させるにはNEGを使います、
入手したGeoTiffを base.tif という名前にして、
gmt grdmath "base.tif" NEG = "sakasama.nc"とすれば、標高が逆転したDEMが作成されます。
SAGAを使ってみる
今回は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.tifとsakasama_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"と書けて、座標系も維持されるので便利です。
今回はここまで。
コメント