gmt grdtrackの使い方

gmt grdtrackとは、任意の地点でのグリッドデータのZ座標を求めるものです。
簡単に言えば、任意の地点のDEMの標高を求めることができます。
gmt project と併せて使うことが多いでしょう。

gmt grdtrack [経度・緯度データ] -G[DEMデータ] -f[座標系] > [出力ファイル名]
※ DEMデータにはgrdcutしたものを指定すること。
※ 地理座標を使う場合はfg。直交座標を使う場合はfc。
※ 出力ファイルは 「経度・緯度・[諸々]・標高」 の順番。

等高線の断面図を描画する

等高線の断面図を描画する

地図に任意の線を描き、その線を通る標高断面図をプロットしてみよう。

今回は富士山周辺の標高断面図を作成してみます。

gmt grdcontourの使い方

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

#!/bin/bash
gmt begin fujisan jpg
    gmt set GMT_VERBOSE = long
    gmt basemap -JM12 -R138.487240/138.954132/35.200372/35.552028 -Bafg -BWSNE
    gmt coast -Df -W0.25 -Lx10/-1+c35+w20+f
    gmt makecpt -Cgeo -T-8000/8000/200 -Z
    gmt grdcut @earth_relief_01s -Gdem.nc -R138.487240/138.954132/35.200372/35.552028 
    gmt grdgradient dem.nc -Ggrad.grd -A45 -Ne0.5
    gmt grdimage dem.nc -Igrad.grd
    gmt grdcontour dem.nc -A500+ap+f6p,0,black -C250 -L0/4000
    gmt plot -W1,red << END
138.487240 35.3625
138.954132 35.3625
END
gmt end

距離を求める

gmt project の使い方 / 震源断面図を描くその3

GMT6のprojectモジュールの使い方の一部を説明しています。

grdtrackで標高を求めるには、座標データが必要なので、projectを用いて0.5km間隔くらいで赤線上の座標を求めます。

#!/bin/bash
gmt project -C138.487240/35.3625 -E138.954132/35.3625 -G0.5 -Q > temp.txt

この結果が、こちらになります。

標高を求める

ここまで来たらgrdtrackを使います。

#!/bin/bash
gmt grdcut @earth_relief_01s -Gdem.nc -R138.487240/138.954132/35.200372/35.552028
gmt grdtrack temp.txt -Gdem.nc -fg > height.txt

この結果が、こちらになります。

1枚の図にする

断面図を描くには、height.txtの1列目と4列目が必要です。
これらを図にすればいいので、

#!/bin/bash
gmt begin fujisan jpg
    gmt set GMT_VERBOSE = long
    gmt basemap -JM12 -R138.487240/138.954132/35.200372/35.552028 -Bafg -BWSNE -Y60
    gmt coast -Df -W0.25 -Lx10/-1+c35+w20+f
    gmt makecpt -Cgeo -T-8000/8000/200 -Z
    gmt grdcut @earth_relief_01s -Gdem.nc -R138.487240/138.954132/35.200372/35.552028 
    gmt grdgradient dem.nc -Ggrad.grd -A45 -Ne0.5
    gmt grdimage dem.nc -Igrad.grd
    gmt grdcontour dem.nc -A500+ap+f6p,0,black -C250 -L0/4000
    gmt plot -W1,red << END
138.487240 35.3625
138.954132 35.3625
END
    gmt basemap -JX12 -R138.487240/138.954132/-100/4000 -Bxagf+l"longitude" -Byagf+l"Height [m]" -BWSne -Y-15
    awk '{print $1" "$4}' height.txt | gmt plot -W1 -Wred
gmt end

とすれば、

となります。