GrADSの出力は、画面に図として表示するのがおもですが、 数値をファイルに書き出すこともできます。 ただし、GrADS (1.9までの群のバージョン)の標準的機能だけでは、 バイナリファイルを書き出すことしかできません。 【2.0群のバージョンではテキストファイルを書き出すこともできますが、 その書式がわたしの意図する目的にはあまり便利ではありません。】 awk, gnuplot, GMTなどで受けつけるテキストファイルを書き出すためには、 GrADSの「スクリプト」機能や、「ユーザー定義関数」機能を使う必要が あるようです。その手続きをいくつか試作したので、順次紹介します。
set gxout fwrite
display (d)コマンド
による出力は、画面表示でなく、grads.fwriteという名前のファイルに
書き出されます。
内容は、実行している機種の内部表現による4バイト浮動小数点のバイナリ数値が
つまっているだけのファイルです。
別のsetコマンドをあらかじめ実行しておくことで、
ファイル名を変えることもできるようですが、わたしはまだためしていません。書き出されたものに、適当な.ctlファイルをつければ、
gradsで読むことができます。
.ctlファイルの書きかたをひとことで言うのはむずかしいですが、
fwriteで書きだされたものは、緯度は南からになっており
(yrevではない)、
レコード長さフィールドはない(sequentialではない)ものです。
【レコード長さフィールドがついているかどうかは、 GrADSのバージョンやオプションによって違う可能性があります。 書き出されたファイルの大きさ(バイト数)を、書き出されるべき数値の個数と 比較してみれば確認できるでしょう。】
series.gsというスクリプトを作りました。
これを、GrADSがスクリプトファイルを読むところ
(環境変数 GASCRP で指定されたディレクトリ) に置けば
run series.gsで使えるはずです。
gradsを起動し、
処理したいファイルをopenしておき、
もし鉛直方向に複数レベルがあれば、その1つに固定したうえで、run series.gs 変数名 経度 緯度 最初のt 最後のt
tとは、
ここでopenしているデータファイルの時刻の通し番号(最初が1)です。
たとえば、xaprecm.ctlから読まれるデータに
月降水量のデータがprecという変数名で収録されていたとします。open xaprecm run series.gs prec 140 35 1 12
gnuplotでplot "prec.lon140.lat35"
なお、
datetoyr.awk
を使うと、
series.gsの出力の第3列を加工し、
1列めに時刻を「年」単位の小数で表わしたものを書き出しますので、
これをgnuplotでプロットすれば、横軸が「年」になります。
aaveser.gs
というスクリプトも作りました。同様にすれば次の形で使えます。run aaveser.gs 変数名 経度西 経度東 緯度南 緯度北 最初のt 最後のt
open xaprecm run aaveser.gs prec 120 150 -20 -10 1 3
prec.120.150.-20.-10というテキストファイルに、
東経120--150度、南緯20--10度の面積平均値の最初の3か月ぶんが
書き出されます。GrADSのユーザー定義関数(udf)として、
wxyvとwxvというものを作りました。
これは、2000年ごろGrADSのFTPサイトに置かれていたUDFの例文 (作者はMike Fiorinoさんらしい)をもとに、わたしが改造したもので、 Fortranで書かれています。 これを、コンパイルして、実行プログラムファイルをある場所に置くとともに、 ユーザー定義関数設定テーブル (udft)に、それを記述する必要があります。 udftを置く場所は、環境変数GAUDFTで指定されます。 これ以上詳しいことが知りたいかたは、わたしにお問い合わせください。 ただし、わたしは2011年4月以後GrADSを実地に使っていませんので、 すぐ応答できるかどうかはわかりません。
【udfという機能は、GrADSのバージョン1.8群や1.9群にあったのですが、
バージョン2.0群ではなくなってしまいました。
udfのプログラムの作りかたはわかりにくく、
わたしにも例文を改造して試行錯誤で動作を確かめることしかできていませんので、
この機能が廃止されるのももっともだとは思いますが、
それに代わる機能が用意されることを希望しております。
なお、バージョン2.0群では、set gxout printとすれば
次のdisplayデータをテキストファイルに書き出すことができますが、
その出力は表示対象の数値が適当な個数ごとに改行されながらならんでいるだけで、
座標値はついていません。】
GMTのxyz2grdが受けつけるxyzファイル、つまり、
各点が1行になっていて、スペース区切りで経度、緯度、変数の値が
ならんでいるテキストファイルを書き出します。
(ただし、実行したときのGrADSの軸の設定によっては、2つの座標は
経度と緯度ではなく、たとえば緯度と気圧かもしれません。)
この前の試作版では名前をwxyzとしたのですが、
GrADSではzは鉛直座標をさすので、valueのvに変えました。
define機能の呼び出しを使います。
書き出したい変数の名前がvならば、define w=wxyv(v)
xyvという名前のテキストファイルができます。
なお、ここで仮にwとした変数にはいる値は、書き出されたファイルの行数(つまり格子点数)
で、'q defval w 1 1' say subwrd(result,3) ' lines written.'
wxyvと基本的に同じですが、
1つの軸の座標値と変数の値の2列だけを書き出します。define w=wxv(v)
xvという名前のテキストファイルができます。vという2次元の変数の東西全経度平均(zonal mean)はdefine vzonalm = ave(v, lon=0, lon=360, -b)
vzonalmは1次元格子点データなので、define w=wxv(vzonalm)
xvに書き出せます
(yvとしたほうがわかりやすいと思いますが、
呼び出された手続きの側では
1次元の座標が緯度かどうかを知らないしくみなので、
こうなっています。)
このxvは、gnuplotでplot "xv"
vの全経度平均を縦軸にしたグラフにしてみる
ことができます。