簡単な解析・可視化
解析・可視化ツールの準備
dcpam5 では入出力するファイルとして Gtool4 NetCDF 規約 に基づいた NetCDF データを扱います.
数値実験の結果を解析・可視化するためには, NetCDF データを取り扱うことのできる解析・可視化ツールが必要です. ここでは, 電脳 Ruby プロジェクト から提供される Gphys を使った可視化の例を紹介します.
可視化ツールのインストール
電脳Ruby謹製品 インストールガイド を参照してください.
GPhys/GGraph による解析と可視化
ここでは, 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 つの分布の描画を行うことができます. mean('lon') は経度方向に単純平均をとることを意味しています.
経度方向, あるいは時間方向には等間隔なデータなので mean('lon') や mean('time') で平均操作ができますが, 緯度方向あるいは高度(圧力)方向には不等間隔なので少し工夫が必要です. 緯度方向に平均操作を行うには緯度格子の重み変数である lat_weight 変数を 読み取り, 重みつき平均をとることにより実現できます.
15: lat_weight = GPhys::IO.open('Temp.nc', 'lat_weight')
16: tempbar = (temp*lat_weight).sum(1)/lat_weight.sum
17: line tempbar.cut('lon'=>0,'sig'=>1)
16 行目で温度に lat_weight 重みをつけた緯度平均をとっています. temp*lat_weight は次元の異なる GPhys オブジェクトの掛け算ですが, GPhys ではきちんと軸を判別して, lat_weight を経度・高度・時間に関して拡張して要素別掛け算をしてくれます. sum(1) は 2 番目の次元に関して和をとることを意味しています. 次元の順番の数え方は 0 番目から始まることに注意してください. lat_weight.sum は重みの総和です. 最後に経度 0, 最下層での温度緯度平均の時間変化を 折れ線グラフで図示しています.
高度平均を計算する場合にはσ座標の重み変数 sig_weight を用いておなじように行うことができます.
このように, GPhys/GGraph を用いると多彩な解析と可視化を実現できます. 処理が長くなってきたら irb でインタラクティブに行うかわりに, エディタを用いてスクリプトファイルを書く方が効率的であり再利用も容易になります. より高度な解析・可視化を行う際には, GPhys チュートリアル を参照してください.
GPhys/gpコマンドによる可視化
ここでは, 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 による解析と可視化が適しているでしょう.
$Id: visualization.rd,v 1.8 2012-03-22 20:49:28 noda Exp $