dcpam5 では入出力するファイルとして Gtool4 NetCDF 規約 に基づいた NetCDF データを扱います.
数値実験の結果を解析・可視化するためには, NetCDF データを取り扱うことのできる解析・可視化ツールが必要です. ここでは, 電脳 Ruby プロジェクト から提供される Gphys を使った可視化の例を紹介します.
電脳Ruby謹製品 インストールガイド を参照してください.
ここでは, Polvani et al. (2004) の傾圧不安定波動実験 で得られたデータを GPhys/GGraph を用いて可視化してみることにします.
まず irb を起動してください.
$ irb
以下のような irb のプロンプトが表示されます.
irb(main):001:0>
このプロンプトに, 以下のようにコマンドを打ちます. 左端の数字は行番号で, 打つ必要はありません.
1: require "numru/ggraph"
2: include NumRu
3: gphys = GPhys::IO.open('Temp.nc', 'Temp')
4: DCL.gropn(1)
5: DCL.sgpset('lcntl', false) ; DCL.uzfact(0.7)
6: GGraph.tone gphys
irb のプロンプトにおいて quit と打つと irb を終了することができます.
ここでは Temp.nc というファイルの中の Temp という変数を読み込み, 図示を行っています.
Temp は経度・緯度・圧力(高度)・時間の 4 次元データですが, 6 行目で図示する際になにも指定しなければ, 後者 2 つの次元に関しては自動的に 1 番目の要素が選択されます. したがってこの図は最下層の初期時刻での緯度経度面の温度を示していることになります
また, 続けて, 下のように時刻を指定することで, 異なる時刻での最下層の温度分布を描くことができます.
7: GGraph.tone gphys.cut('time'=>10)
終了させるには
8: DCL.grcls 9: quit
としましょう.
絵を描くだけでなく, 解析を行うこともできます. 例として温度と圧力座標から温位を計算し, 図示してみましょう.
1: require "numru/ggraph"
2: include NumRu
3: temp = GPhys::IO.open('Temp.nc', 'Temp')
4: sig = GPhys::IO.open('Temp.nc', 'sig')
5: ps = GPhys::IO.open('Ps.nc', 'Ps')
6: kappa = 2.0/7.0
7: theta = temp*(ps/1e5)**kappa/(sig**kappa)
8: DCL.gropn(1)
9: DCL.sgpset('lcntl', false) ; DCL.uzfact(0.7)
10: GGraph.tone theta.cut('time'=>0,'lon'=>0)
11: GGraph.contour temp.cut('time'=>0,'lon'=>0),false
12: GGraph.color_bar
第 7 行目で温位の計算を温度と圧力から行っています. 図示させているのは経度 0 度の子午面断面です. 温度を等値線で, 温位を色で塗分け表示しています.
1 日目と 10 日目の温位の経度平均の違いを見たければ, 続けて
13: GGraph.tone theta.cut('time'=>0).mean('lon')
14: GGraph.contour theta.cut('time'=>10).mean('lon'),false
と入力することで同時に 2 つの分布の描画を行うことができます.
このように, GPhys/GGraph を用いると多彩な解析と可視化を実現できます. 処理が長くなってきたら irb でインタラクティブに行うかわりに, エディタを用いてスクリプトファイルを書く方が効率的であり再利用も容易になります. より高度な解析・可視化を行う際には, GPhys チュートリアル を参照してください.
ここでは, Polvani et al. (2004) の傾圧不安定波動実験 で得られたデータを GPhys 付属の gp コマンドを用いて可視化してみることにします. 温度のデータを読み取り図示するには,
% gpview Temp.nc@Temp
と入力します. これは Temp.nc というファイルの中の Temp という変数を読み込み, 図示せよというコマンドです.
Temp は経度・緯度・圧力(高度)・時間の 4 次元データですが, なにも指定しなければ後者 2 つの次元に関しては自動的に 1 番目の要素が選択されます. したがってこの図は最下層の初期時刻での緯度経度面の温度を示していることになります. 断面を変えて, 時刻 10 日後の最下層の図にしたければ
% gpview Temp.nc@Temp,time=10,sig=1
と, カンマで区切って断面を指定することができます.
時間の代わりに子午面断面などを描きたければ
% gpview Temp.nc@Temp,lat=0,time=10
経度平均操作もできます. ただし緯度・高度方向は不等間隔なので, 単純平均操作は間違ったものになってしまうことに注意してください.
% gpview Temp.nc@Temp,time=10 --mean lon
アニメーションも簡単に見ることができます. --Gaw オプションをつければページ送りが自動的に行われます. 以下は最下層の温度分布の時間発展を見るコマンドです.
% gpview Temp.nc@Temp,sig=1,lat=10:80 --anim time % gpview Temp.nc@Temp,sig=1,lat=10:80 --anim time --Gaw --wsn 4
lat=10:80 は緯度範囲を 10 度から 80 度までに限定するオプションです. --wsn 4 は DCL の装置番号を指定するオプションです(4 は GTK).
座標を 3 つ指定して 1 次元データにすると折れ線グラフを描けます.
% gpview Temp.nc@Temp,lon=0,lat=0,time=0 --exch
--exch オプションは縦軸と横軸を入れ換える操作を指示するものです.
gpview には他にもいろいろなオプションがあります. gpview --help とするとオプションと使い方の例が表示されます. gp コマンドシリーズには他にもいろいろなものが用意されています. 主なものは以下の通りです.
gp コマンドシリーズは 1 行入力ですぐに結果を表示できるのが特徴です. そのためクイックルックや計算の途中でのデータチェック等に便利です. しかしながら複数の変数を組み合わせた解析や可視化はできません. 本格的なデータ解析と可視化を行うには, GPhys/GGraph による解析と可視化が適しているでしょう.
