gmt grdgdalの使い方

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

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

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コマンド

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

gdaldem：-Adem+m【hillshade, color-relief, slope, TRI, TPI roughness】
DEMであれこれ遊べます。陰影図とかも作れます。DEMのカラーパレットの変更もこれ。

gdaldem — GDAL documentation
gdal.org
gdalinfo：-Ainfo
ラスタやベクタファイルの情報を表示できます。

gdalinfo — GDAL documentation
gdal.org
gdal_grid：-Agrid
散布状に存在している点データを格子状のデータに最近傍等で格子状データに補完します。

gdal_grid — GDAL documentation
gdal.org
gdal_rasterize：-Arasterize
ベクタファイルをラスタライズします。（そのまんまですが…）

gdal_grid — GDAL documentation
gdal.org
gdal_translate：-Atranslate
ラスタファイルを色々変換できます。フォーマットを変換したり、座標系を変えたりもできます。

gdal_translate — GDAL documentation
gdal.org
gdalwarp：-Awarp
座標変換や地図投影法の変換ができます。

gdalwarp — GDAL documentation
gdal.org

地図投影法を変換する

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で地図投影法の変換を行ってみます。
-Aにwarpとだけ入れて残りのオプションは-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でもうまく行っているような気がします。

座標を確認する（gmt grdinfo）

gmt info・gmt grdinfoの使い方

GMTのgrdinfo・grdinfoモジュールについて説明しています。

gmt grdinfoを使えばもっと便利にできます。

変換後の画像にも座標の情報が記されているので、gmt grdinfoで確認しておきましょう。

gmt grdinfo nagata.tiff

としてみると、

nagata.tif 139.732249882 139.759379924 35.6694925022 35.686881192 0 255 8.37605481785e-06 8.37605481785e-06 3239 2076 1 1

と表示されるので、これを配列に格納すれば便利ですね。

#!/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"
    range=(`gmt grdinfo nagata.tif  -C`)
    gmt grdimage nagata.tif  -JM12 -Bafg -BWSNE -R${range[1]}/${range[2]}/${range[3]}/${range[4]}
gmt end

とすればすっきりですね！

陰影図を作成する

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

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

#!/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

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

・・・微妙ですね。