[ GFD Dennou Ruby Home | GPhys Home | Tutorial top ]

GPhys, GGraphチュートリアル (その7)

GPhysを使ったデータ解析*

Ruby と GPhys を使いこなせば,地球流体の研究で行う様々なデータ解析ができます. 本章はそのための入門編です.この先,GPhysレシピ を参考にさらに経験を積めばなんでも (^_^) できるようになります. *1

目次

平均などの統計操作

ある座標軸に沿って平均するという操作は,これまでも出てきました. ここではより包括的に解説します. mean は,格子点に関するアンサンブル平均です (格子間隔にもとづく重み等はかけません). その仕様は次のようになってます.

GPhys#mean(*dims)

ここでは Ruby 流の記述をしました.これにも慣れましょう:

mean と同じ仕様のメソッドには次があります:

いずれも NArray, NArrayMiss のメソッドです.データに欠損があるときは NArrayMiss が使われます. NArrayMiss は,欠損は飛ばして処理します(平均するときの母数にもカウントしません). GPhys は処理の中身はこれら委ねます.

では,サンプルプログラムです.

実行してみましょう.

$ ruby sample_mean_etc.rb 
gp :
<NumRu::GPhys grid=<3D grid <axis pos=<'lon' shape=[144]  subset of a NumRu::VArrayNetCDF>>
        <axis pos=<'level' shape=[8]  subset of a NumRu::VArrayNetCDF>>
        <axis pos=<'time' shape=[31]  subset of a NumRu::VArrayNetCDF>>>
   data=<'air' shape=[144, 8, 31]  subset of a NumRu::VArrayNetCDF>>
pressure levels:
NArray.sfloat(8): 
[ 925.0, 700.0, 500.0, 300.0, 200.0, 100.0, 50.0, 20.0 ]

gp.mean(0) :
<NumRu::GPhys grid=<2D grid <axis pos=<'level' shape=[8]  subset of a NumRu::VArrayNetCDF>>
        <axis pos=<'time' shape=[31]  subset of a NumRu::VArrayNetCDF>>>
   data=<'air' sfloat[8, 31] val=[276.323852539062,266.1376953125,250.100402832031,226.781204223633,...]>>

gp.mean('lon','time') :
<NumRu::GPhys grid=<1D grid <axis pos=<'level' shape=[8]  subset of a NumRu::VArrayNetCDF>>>
   data=<'air' sfloat[8] val=[274.950927734375,264.612152099609,248.710037231445,224.867828369141,...]>>

gp.mean :
NumRu::UNumeric
234.546216957885 degK

mean,max,min,stddev values:
NArrayMiss.sfloat(8):
[ 274.951, 264.612, 248.71, 224.868, 217.77, 214.26, 213.517, 217.683 ]
NArrayMiss.sfloat(8):
[ 287.55, 278.02, 261.12, 235.8, 236.57, 226.35, 226.35, 235.0 ]
NArrayMiss.sfloat(8):
[ 253.87, 245.42, 230.77, 215.32, 205.35, 202.62, 200.57, 200.32 ]
NArrayMiss.sfloat(8):
[ 6.73178, 6.87236, 6.20152, 4.12792, 6.29422, 4.00543, 4.48009, 6.30728 ]

ソースの最初は require "numru/gphys" となってます.可視化を行わない場合,numru/ggraph でなく numru/gphysrequire することで起動が少し早くなります*2

プログラムでは,最初に gp.cut.... の行でサブセットをとります. GPhys のサブセット切り出し([]cutによる)は, もとのデータのサブセットへのマッピング法を格納したオブジェクトを作るだけなので, この時点ではデータは読み込みません(だからもしもサブセットに書き込めば参照元に書き込まれます). このことは,出力の最初のほうに subset of a NumRu::VArrayNetCDF と現れることから伺えます. データ読み込みは gp.mean(0) などの演算を行う際にはじめて行われます. その際は必要なサブセットだけが読み込みまれます. ですから,サブセット切り出しは演算の前に行うのがよいです.

結果をみると,mean の戻り値は引数ありの場合 GPhys, なしの場合 UNumeric であることがわかります(省略せずに書くと NumRu::GPhys と NumRu::UNumeric です).gp.mean(0) の結果の表示では, data=<'air' sfloat[7] val=[274.950927734375,264.612152099609,... となってることに注目してください.演算の結果はメモリー上に NArray または NArrayMiss として置かれます.つまり,同じ GPhys オブジェクトでも, 最初はデータはファイル中にありますが,演算の結果できたものは全体が計算機のメモリー上に保持されるのです. このように違っていても同じ GPhys なので同じように描画できるというのは, これまでの章で体験しました(こういうことが簡単に実現できるのがオブジェクト指向言語の利点です). なお,演算結果をファイルに書き出す方法は後で示します.

プログラムの最後では,緯度40度の各気圧面から得られたいくつかの統計量を表示します. 母集団は経度および時間の2次元断面にかかる 144×31=4464 個のデータです. 入力データは欠損があり得る形で NetCDF ファイルに格納されているので 結果は NArrayMiss クラスですが,実際には欠損はありません (はじめる前にの末尾を参照).

事例:半旬平均で見るブロッキング高気圧

次のプログラムは2012年1月後半の500hPa等圧面高度の半旬(5日)平均値(コンター)と, その気候値からの偏差(色)を表示します. 最初の図(1月16-20日)にみられる日付変更線付近の高緯度の高気圧偏差はブロッキング高気圧です. タイトルに平均期間が表示されるよう工夫してあることに注意してください.

結果の一枚: img/sample_pentad_s.png(クリックでフルサイズ表示)

演習問題

  1. 上記の sample_pentad.rb を編集して, 異なる期間や異なる時間間隔で作図してみましょう. また,半旬平均をとらず日平均データを表示してみましょう.

ファイルへの出力

データ処理した結果をファイルにとっておきたいことは多いでしょう. GPhys では次のように簡単に NetCDF ファイルに書き出せます.

実行と結果確認:

$ ruby sample_output.rb 
[horinout@localhost trunk]$ ncdump -h air.2012-01.monmean.nc 
netcdf air.2012-01.monmean {
dimensions:
        lon = 144 ;
        lat = 73 ;
        level = 17 ;
variables:
        float lon(lon) ;
...中略...
        float air(level, lat, lon) ;
...後略...

出力のポイントは3つ:

です.GPhys::IO.write は,座標軸を含む GPhys オブジェクトを構成する要素をすべてファイルに書き出すので,これだけで自己記述性のあるデータが書き出されます.

なお,ファイル中の名前は重複しないようにしなければなりません. よって,一つのファイルに同じ名前の複数の変数(例えば "air" の経度平均と時間平均) を出力したければ少なくともいずれかは GPhys"name= メソッドを使って改名する必要があります.

地球流体力学・気象学的なデータ解析

GPhys には他にもいろいろな数値計算手法, さらに気象学・地球流体力学的なデータ解析手法を提供するライブラリがあります. 数値微分,フーリエ解析,球面幾何,気象学的な諸量(温位,渦位等)の計算等々です. それらについては, GPhysのマニュアル のほか, GPhysレシピ でも紹介する予定です. 以下では.GPhys付属のライブラリ利用例を示し,自分でライブラリを作る例を示します.

GPhys が提供するライブラリ.事例:渦位計算

次の例は,渦位(ポテンシャル渦度; potential vorticity; PV)を計算して図示します.

最初の行は,これまでのソースにはない require "numru/ganalysis" となっています. こうすると GPhys において応用的なライブラリを集めた GAnalysis というモジュールが利用できるようになります. GAnalysis は GPhys パッケージの一部ですが, このように陽に指定しない限り読み込まれないようになっています. これにより,GAnalysis が必要ない場合の GPhys の立ち上げが速くなってます.

上のプログラムでは渦位は pv = GAnalysis::Met.pv_on_theta(u, v, theta, theta_levs) によって計算されます.これは GAnalysis 中の Met というモジュール中の pv_on_theta というメソッドで,気圧座標の入力データを元に等温位面(一般に複数個)上での渦位を計算します (GAnalysis::Metのオンラインマニュアル参照). 4つの引数は順番に,東西風,南北風,温位(以上同じ格子点上で定義された GPhys オブジェクト),渦位を求めたい温位面(数値のArray)です.

プログラム冒頭に書かれているように,4つのコマンドラインオプション(引数)があります. 最初の引数では温位面を指定します.以下はデフォルトの 320 K での結果例です. 上述の sample_pentad.rb で表示したブロッキング高気圧の形成過程を示します. 低緯度の低渦位の空気塊が高緯度に伸び切離したことがわかります.

実行結果より:

img/sample_pv_0114_s.png img/sample_pv_0115_s.png img/sample_pv_0116_s.png img/sample_pv_0117_s.png img/sample_pv_0118_s.png img/sample_pv_0119_s.png img/sample_pv_0120_s.png

演習問題

  1. 上記の sample_pv.rb にコマンドライン引数を与えて, 温位 700 K での渦位の時間発展をみてみましょう. また,表示する期間を変えたり,アニメの時間間隔を変えてみましょう (0を指定するとマウスクリックを待ってから次のページを表示するようになります).
  2. 上記の sample_pv.rb をもとに,1ページに2枚画像を 表示できるようにし,渦位と等圧面高度を表示するようにしましょう (等圧面高度は sample_pentad.rb と 同じものを時間平均せずに表示しましょう).

自分でライブラリを作る.事例:Takaya&Nakamura (2001)の波活動度フラックス

さて,付属のライブラリがどんなに充実してもすべての研究者の需要はみたせません. データ解析手法を自分で実装しライブラリ化して使いまわせるようになりましょう. ライブラリ化といっても Ruby では簡単で,独立したファイルにするだけです. そのファイルを require すれば利用できます. (require については,一般的な Ruby の参考書等を参照してください). 以下に作成例をひとつ示します.サンプルプログラム付きのライブラリになります.

とりあげるのは Takaya&Nakamura (2001)の波活動度フラックスです. それは次の式で表されます (Takaya and Nakamura, 2001, J.Atmos.Sci, 608-627, Eq.(32)).

img/TN2001eq32.png

ここで,U = [U,V,0] は「基本場」の東西,南北流,ψ' は準地衡流線関数の偏差,p は気圧を 1000 hPa で規格化したものです. x,yは東西,南北で球面の効果を無視しています (球面版の実装も難しくありませんが,ここでは簡単のため x,y 版にします). Cu は波の位相速度に関連する量で停滞波ではゼロになるので,ここでは Cu M の項は無視します(いつでも無視できるわけではないので適用限界があります). なお,簡単のため,水平成分だけ計算することにします.

この波活動度フラックスの計算を実装してみます. (慣れればこの程度の複雑さのものはすぐ作れます.しかも GPhys メソッドとして作るため入力に関しかなり汎用になります.) では,ちょっと長いですがソースをざっと眺めてみましょう. ライブラリの主な(実質的に唯一の)部分はメソッド tn2001fluxCSH です. 中をみると,数式との対応が見て取れるのではないでしょうか. (なお,上記の事情を勘案して背景風の取り扱いについて二通りの計算が用意してあります.)

メソッド tn2001fluxCSH

module NumRu
  class GPhys
  end
end

にくるまれてるので,GPhys のメソッドとして定義されていることがわかります. 微分は差分で近似しました.経度,緯度座標はあらかじめ x = a cosφ0 λ, y = aφ に変換しておきます(aは地球半径,λは経度,φは緯度,φ0は基準となる緯度です). 使った差分法は2次精度で不等間隔格子に対応します(一階微分は等間隔格子なら中央差分になります). 仕様についてはマニュアル(Reference manuals of the Derivative library) を参照してください.

ライブラリ定義の下の部分()は

if __FILE__ == $0
end

で囲われています.このファイルを他のファイルから require するとこの中は実行されず,このファイルを直接 ruby コマンドの引数とした場合のみ実行されます. ライブラリの利用法のサンプル兼テストプログラムになっている次第です. このやり方は,Ruby でライブラリを作る時の常道の一つです *3

では実行してみましょう.

% ruby tn2001flux.rb

結果の一枚: img/tn2001flux_s.png(クリックでフルサイズ表示)

参考:背景風(最終ページ): img/tn2001flux_climUV_s.png(クリックでフルサイズ表示)

ここでは背景場は気候値にしました. 本計算には停滞性の仮定があるため,先程より長い10日間の平均をとって, 気候値からのずれに対し波活動度フラックスを計算しました. 左右の図の違いは矢印で表した波活動度フラックスの計算法のみです. 左は Cu = 0 とした式(32)そのままですが, 右は背景南北風がゼロと仮定した場合に相当します. 左右が大きく異なるところは,|V| が |U| に比して小さいと言えないわけですが, 一般に |V| は小さいので,背景風自体が小さい可能性があります(上図のブロッキング高気圧あたりがそうです). そうなると,背景風に比して位相伝播を無視するという停滞性の仮定はよくないかもしれないので注意が必要です(なお,今回無視した CuM を計算するには波の位相速度を特定しないとなりません). ちなみに,今回のように CuM を無視すると,波の位相伝播から相対的にみた波活動度のフラックスになります.その意味で,停滞性の波でなくても意味がある計算ではあります.

なお,本例に対し論文の著者の高谷康太郎さんより次のようなコメントをいただきました.

”アラスカ付近のブロッキングを対象とされていますが、例としては一番解釈の難しいケースを取り上げられたな、という感じはします。御指摘の通り、アラスカ付近は気候値でも背景場が歪んでおり背景風速の弱い領域です。...(中略)...チュートリアルとしては、むしろ、ジェット中のブロッキングであるヨーロッパ付近やユーラシア大陸上のものの方が良かったかも知れません。もちろん、ブロッキング高気圧そのものは非線形現象ですが、ブロッキングへの上流からのロスビー波の入射、またはブロッキングから下流への射出、などは線形性の強い現象であり、フラックスの例としてはより良いかと思います(TN2001またはTakaya and Nakamura 2005など)。”


*1「なんでも」というのは言いすぎかもしれませんが,GrADSなどと違い Ruby は本格的なプログラミング言語なので自由度は非常に高いです. 計算速度がネックになれば C で拡張することで補えますし.
*2numru/ggraphrequire するとその中で numru/gphysrequire されます
*3もっと厳格に「テスト駆動開発」もできますが,ここでは一番簡単なやり方として紹介します.いずれにしてもサンプルは欲しいですし.


davis Group / GFD Dennou Staff dcstaff@gfd-dennou.org