研究していて見つけたちょっと便利そうなことなど書いていこうと思います。
GMTを使うとpsファイルとしていろいろときれいな図を作ることができますが、波形を描こうとするときはちょっと不便です。というのも-Rオプションで範囲指定をしなければならないので、あらかじめ最大値を知らなければならないからです。私はawkを使って最大値を読み取っています。
1 0.328476E-01
2 0.328479E-01
3 0.328484E-01
・ ・ ・
上のような波形ファイルの1401行目から3400行目を読み取って波形を描くならば、
set file_in = hoge.dat
set start = 1401
set end = 3400
set max = ` awk 'BEGIN {max=0} \
$1 >= '"$start"' && $1 <= '"$end"' && \
$2*$2 > max*max {max=$2} \
END {if (max > 0) print max; \
else print -1*max}' $file_in `
psxy $file_in -JX20/8 -R${start}/${end}/-${max}/${max} （以下略）
とすると、きちんと枠いっぱいに収まった波形を描くことができます。同じファイルを２回読み込むことになるため少々時間はかかりますが、いちいち他のプログラムなどで最大値を読み取る必要がいらないので個人的には重宝しています。Cシェルの配列を使えば複数のファイルをがしがし描き出すこともできますし。(2006/08/27)
$ cat hoge1.ps hoge2.ps hoge3.ps > hogehoge.ps $ ps2pdf hogehoge.ps
catコマンドでpsファイルも結合できることを知ってちょっとびっくり。ただしWINDOWSのghostviewでは複数ページの表示はできませんでした。（2006/09/08）
「GMT 波形」とか「Cシェル 配列」などという検索で当ページにたどり着く方もいらっしゃるようなのですが、それについてはほとんど触れていなくて申し訳ないので、上記の２つのネタも含めて波形を描いてみましょう。
#!/bin/csh -f
set xshift = 9.5
set dir = nj06
set stn = ( HRS003 HRS008 HRSD06 OKY005 OKY008 OKYD02 OKYD03 OKYD05 OKYD13 OKYD14 SMN002 SMND10 TTR005 TTR006 TTR008 TTRD04 )
@ istn = 1
@ j = 1
while ( $j <= 4 )
set file_out = 4stn${j}.ps
echo "-2.0 3 24 0 0 CM invtt2_ini1multi4skp4 $dir " | pstext -JX21/5 -R0/21/0/5 -N -K -X12.5 -Y25.5 -P >! $file_out
echo "-9.5 2 20 0 0 LM EW-component" | pstext -J -R -N -O -K >> $file_out
echo " 0.0 2 20 0 0 LM NS-component" | pstext -J -R -N -O -K >> $file_out
@ i = 1
while ( $i <= 4 )
set file_in = ${dir}/${stn[$istn]}.ew
set max = ` awk 'BEGIN {max=0} \
$2*$2 > max*max {max=$2} \
END {if (max > 0) printf("%.2f", 100*max); \
else printf("%.2f",-100*max)}' $file_in `
echo "Maximum velocity in $file_in is $max"
set yrang = ` echo "scale = 2; $max*2" | bc`
psbasemap -JX8/4 -R0/2048/-${max}/${max} -B/a${max}W -O -K -X-${xshift} -Y-4 >> $file_out
echo "20.48 5.5 18 0 0 CM $stn[$istn]" | pstext -J -R0/20.48/-4/4 -N -O -K >> $file_out
awk '{print $1,$2*100}' $file_in | psxy -JX8/8 -R0/2048/-${yrang}/${yrang} -W3 -O -K -Y-2 >> $file_out
awk '{print $1,$3*100}' $file_in | psxy -J -R -W3/255/0/0 -N -O -K >> $file_out
set file_in = ${dir}/${stn[$istn]}.ns
set max = ` awk 'BEGIN {max=0} \
$2*$2 > max*max {max=$2} \
END {if (max > 0) printf("%.2f", 100*max); \
else printf("%.2f",-100*max)}' $file_in `
echo "Maximum velocity in $file_in is $max"
set yrang = ` echo "scale = 2; $max*2" | bc`
psbasemap -JX8/4 -R0/2048/-${max}/${max} -B/a${max}W -O -K -X$xshift -Y2 >> $file_out
awk '{print $1,$2*100}' $file_in | psxy -JX8/8 -R0/2048/-${yrang}/${yrang} -W3 -O -K -Y-2 >> $file_out
awk '{print $1,$3*100}' $file_in | psxy -J -R -W3/255/0/0 -O -K >> $file_out
@ i++
@ istn++
end
psbasemap -JX8/1 -R0/20.48/0/1 -Ba5/S -O -K -X-${xshift} -Y1 >> $file_out
psbasemap -J -R -Ba5/S -O -X$xshift >> $file_out
\rm .gmt*
@ j++
end
cat 4stn?.ps > 4stn.ps
echo "Now converting PS file to PDF file"
ps2pdf 4stn.ps
これより以下のような波形が描けます。
解説はそのうち書きます。たぶん。（2006/12/12）
ps auxすれば、誰かがバックグラウンドで流しているジョブも含めて表示してくれますが、どうでもいい種々の仕事まで出てくるので一体今アクティブなのはどれなのかが非常に分かりにくいものです。私は次のコマンドでそんな余計なものをそぎ落としてます。
$ ps aux | grep -e "R " -e "R+" | grep -v -e grep -e ps
これを.tcshrcファイルでdareと名づけて使っています。時々最後がRで終わるオプションを持つ空ジョブなどが表示されてしまいますが、ご愛嬌ってことで。topコマンド使えばいいじゃん、って言われるとまあ。 (2007/05/02)
<追記>StatusがR+の場合を検出できないことが判明したので改良しました。 (2007/05/12)