2012年 8月 14日 掲載
今回は，その53〜56であつかった grdgradientコマンドの別の使い方を紹介する．
grdgradientコマンドでは，-Sオプションで指定したファイルに，傾斜の大きさが書きこまれる． このとき，傾斜の方向をさがす -Dオプションも同時に指定する必要がある．
ではさっそくやってみよう． 元の grdファイルには，その53 でも使った人工的な地形 crater.grd を使う（地形を作成する awkスクリプトは crater.awk）．
$ cat crater.awk
BEGIN{
for(x = -5; x <= 5; x += 0.03125){
for(y = -5; y <= 5; y += 0.03125){
r = sqrt(x^2 + y^2)
if(r < 1.5){
z = 0
}else if(r < 2.5){
z = r - 1.5
}else if(r < 3.5){
z = 3.5 - r
}else{
z = 0
}
print x, y, z
}
}
}
$ awk -f crater.awk > crater.txt
$ xyz2grd crater.txt -Gcrater.grd -R-5/5/-5/5 -I0.03125
$ grdgradient crater.grd -Ggradient.grd -A0
$ makecpt -Ccopper -Z > cp.cpt
$ grdimage crater.grd -Jx1 -Igradient.grd -Ccp.cpt -Ba2 > crater.eps
このクレーターの斜面の傾斜を計算するにはこうする．
$ grdgradient crater.grd -Smag.grd -D
できあがった mag.grd を図化してみよう．
$ makecpt -Cpolar -Z > p.cpt $ grdimage mag.grd -Jx1 -Cp.cpt -Ba2 > mag.eps
ちゃんと斜面の傾斜が 1（赤）になっている．
傾斜の方向は，-Gオプションで指定されたファイルに書きこまれる．
$ grdgradient crater.grd -Smag.grd -Garg.grd -D $ makecpt -Ccyclic -Z > c.cpt $ grdimage arg.grd -Jx1 -Cc.cpt -Ba2 -K > arg.eps $ psscale -D5/-1/10/1h -Cc.cpt -Ba90 -O >> arg.eps
傾斜の方向は，上っていく方向が方位（北が 0度，時計回り正）で示される．
[ KDOホームにもどる] [ いちからはじめる GMT インデックスにもどる ]
(C) 2012 KDO