gmt grdgdalの使い方

grdgdal — GMT 6.3.0 documentation

gmt grdgdalとはGMT上でGDALのコマンドを実行するためのモジュールです。
実行には勿論GDALがインストールされている必要があります。

GDALのGMTラッパーと捉えるのがいいかもしれません。
普通にGDALいじった方が便利な気もしますが。。どうなんでしょうか。

GDAL — GDAL documentation

gmt grdgdal 【DEM等のラスタイメージもしくはベクターファイル】 -A【実行したいGDALコマンド】[+m【gdaldemのオプション】][+c【カラーパレット】] -G【出力ファイル名】 -F【gdalコマンドのオプション】 -M【+r/w】
gdaldemを使用する場合はオプションを選択する。オプションは【hillshade, color-relief, slope, TRI, TPI roughness】の中から選ぶこと
-M+rとするとGDAL経由で読み込む。-M+wとするとGDAL経由で出力する。-M+r+wでGDALで読み書き。

使用できるGDALコマンド

現在、使用できるGDALのコマンドは以下の通りです。

地図投影法を変換する

Google MapはEPSG:3857でXYZ tilesを提供していますが、GMTではEPSG:4326あたりにしてあげないとgrdimageがうまく行きません。

まずは適当にGoogle Mapsや地理院地図でGeoTiffを持ってきます。

地図出典:国土地理院

GDALで変換する

これはGDALでやるほうが簡単な気がします。

#!/bin/bash
gdalwarp -s_srs EPSG:3857 -t_srs EPSG:4326 kokkai.tif nagata.tif
gmt begin kokkai png
    gmt grdimage nagata.tif  -JM12 -Bafg -BWSNE -R139.7345/139.7568/35.6673/35.6847
gmt end

で変換・描画できます。

GMTで変換する

grdgdalで地図投影法の変換を行ってみます。
-Awarpとだけ入れて残りのオプションは-Fの中に入れればいいみたいです。

読み込みオプションは-M+rとするとうまく行きました。理由は分かりません。

#!/bin/bash
gmt begin kokkai png
    gmt set GMT_VERBOSE = long
    gmt grdgdal "./kokkai.tif" -M+r -Awarp -Gnagata.tif -F"-s_srs EPSG:3857 -t_srs EPSG:4326"
    gmt grdimage nagata.tif  -JM12 -Bafg -BWSNE -R139.7345/139.7568/35.6673/35.6847
gmt end

とすると、

地図出典:国土地理院

となったのでgrdgdalでもうまく行っているような気がします。
ただ、-Rを指定しないとエラーになるのが気になります。

陰影図を作成する

需要があるか分かりませんが、ここのDEMを取得して陰影図を作成してみましょう。

GDALで陰影図を作成するにはgdaldemhillshadeを用いるので、

#!/bin/bash
gmt grdcut @earth_relief_01s -R139.7345/139.7568/35.6673/35.6847 -Gkokkai.nc
gmt grdgdal -Adem+mhillshade kokkai.nc -Gkokkai2.nc -M+r -F"-of GMT -b 1 -z 1.0 -s 1.0 -az 315.0 -alt 45.0"
gmt begin kokkai2 png
    gmt grdimage kokkai2.nc -JM12 -R139.7345/139.7568/35.6673/35.6847 -Bafg -BWSNE -Cgeo
gmt end

とすれば陰影図を作成できます。

・・・微妙ですね。


そういえばGMT6の日本語解説ブログを作成されている方がいました!分かりやすいです。

私が未だに習得していないgmt moviegmt convert等の使い方も掲載されているので是非見てみてください。

GMTで描こう Part 5: movie | ゆるい趣味ノート

コメント

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