GrADS, gnuplot, GMTを使った時系列のグラフ作図の例として、
毎月の降水量
(アメリカ合衆国NOAA NCEP CPCのXieたちが編集したCMAPと呼ばれるもの)を
赤道付近のおもに海上のいくつかの領域で面積平均した値の時系列をグラフにする。
以下にあげるプログラムは、CNSの
/pub/sfc/earthobs/example/sergraph/にある。
入力データはGrADSで読める形(バイナリ)にしてある
(別記事参照)。
面積平均せず、1つの格子点の値を拾ってグラフにするだけならば、GrADSの displayコマンドだけでできる。 複数の格子点の値を同じ軸上に重ねる例を dispser0.ge として用意した。
グラフを重ねる場合、set vrangeでy軸の値の範囲を設定しておく必要がある。 そうしないと、データごとに別々にスケールされるので、 どのグラフがどの軸に対応するのか読めなくなってしまう。
実行方法は、grads中で exec dispser0.ge である。 これは、赤道上、東経50度から西経70度まで20度おきの各点 (ただしここでは海上の格子点に限りたかったので、 東経110度は陸にぶつかるため省略)の値を重ねて表示する。
これだけたくさん重ねても、見分けることができないので、 実用的には、時系列4個以下くらいにとどめるべきである。
面積平均は、GrADSの関数 aave を使ってできる。 しかし、残念ながら、display aave(...)でグラフを作図することはできない ようだ。aaveを計算するためには、x, y座標がどちらも幅をもって与えられて いなければならないが、この状態でaaveの結果をdisplayすると、グラフに ならず、テキスト画面に数字で出力される。 グラフにするには、define機能を使ってaaveの結果をいったん変数に入れ、 set lonとset latで経度と緯度にひとつずつの値(なんでもよい)を与えた うえで、defineしておいた変数をdisplayすればよい。 このような方法による例を dispser.ge として用意した。 実行方法は、grads中で exec dispser.ge である。
なお、降水量は時間変化のはげしい変数であり、
月単位でみても時系列はなめらかではない。
長期の変化をみるには、数値的に平滑化してみるのがよい。
その方法のひとつとして、移動平均がある。
GrADSでは、たとえばprecma = ave(prec, t-5, t+5)
とすれば、precの11項(例のデータでは11か月)移動平均をprecmaとして
使うことができる。
もし1年周期現象を消したいとすれば、12か月移動平均がよいのだが、
偶数項の移動平均では、代表時刻が中央からずれてしまうという問題がある。
そこで、例文dispsers.geでは、
13項を使い、最初と最後の項の重みが半分になるような移動平均を
GrADSのsum関数と演算式を使って計算し、displayで表示するようにした。
GrADSでは、vpageという機能を使って、1つの画面に複数のグラフを配置する ことができる。vpageでは使う画面の範囲を指定するが、インチ単位で 11×8.5 (縦長の場合は8.5×11)のうちでなければならない。また、グラフは vpageで設定されたうちで余裕をもって配置されるので、画面を有効に使うため には、vpageで設定する領域を相互に重なるように与える必要があるようだ。 (この点はわたしのGrADSの画面サイズ指定に関する理解不足かもしれない。) 4つのグラフを上下にならべる例を dispser4.ge として用意した。 実行方法は、grads中で exec dispser4.ge である。 結果をwi dispser4.pngでPNG画像ファイルにしたものがdispser4.pngである。 (横軸の1978は、1978年の初めを意味する。データの初めは1979年1月である。)
gnuplotやGMTでグラフを作るためには、時系列データをテキストファイルにして おく必要がある。 面積平均値の時系列をテキストファイルとして書き出すためのGrADSスクリプト aaveser.gs が作ってある。 これを利用して、2の例と同じ領域の平均値を書き出す手続きがmakeser.ge、 実行方法は、grads中で exec makeser.ge である。 この結果、prec.100.120.-10.10のような名前のファイルができる。 内容はスペース区切りで、1列めは時刻の通し番号、2列めが(この場合降水量の) 値、3列めがGrADSの表現による日時となっている。
このままだと時間めもりがわかりにくいので、3列めの日付を加工し、 出力ファイルの1列めとして、「年」を単位とする(小数の)日時を 書き出すプログラムがdatetoyr.awk、 これを一連のファイルについて実行する手続きがdodatetoyr.shである。 実行方法は、シェルプロンプトから sh dodatetoyr.sh である。 この結果、PREC.100E.120E.10S.10Nのような名前のファイルができる。
3の処理結果をgnuplotで作図する手続きがtimeser.gpltである。 実行方法は、gnuplot中で load "timeser.gplt" である。
ここでは、清書は次のGMTでやることを想定し、 データの内容を確認するために、 gnuplotで全部まとめて作図してしまったのだが、 同じ軸上に重ねることしかできないので、 データを見て理解できるグラフとはいえない。 実用的には、 点グラフと線グラフを使い分け、それぞれ4つ、 合わせて8つの時系列が限度だろう。
この画面を次の手続きでPNG画像ファイルにしたのがtimeser.gplt.pngである。
set term png set output "timeser.gplt.png" replot
3の処理結果をGMTで作図するシェルスクリプトtimeser.cshである。 GMTでは、軸の位置を自由にずらすことができるので、ここでは、紙上で2 cm (降水量10 mm/日に相当)ずつずらして作図してみた。 結果のPostScriptファイルがtimeser.ps、 それをImageMagickのconvertコマンドで
convert timeser.ps timeser.gmt.pngとしてPNG画像に変換したものが、timeser.gmt.pngである。