GFD Dennou Club / Dennou Ruby / Tutorials

Title: 気象解析ライブラリの開発方法

塚原大輔, 堀之内武

最終改訂: 2005 年 03 月 23 日


1 目次


2 はじめに

本チュートリアルでは GPhys::EP_Flux 内部で用いられている微分演算モジュールを題材に NArray および GPhys を利用した Ruby による気象解析モジュールの開発の仕方を 順を追って紹介していきます.

各節の末尾に関連する小技(TIPS)へのリンクを提示します. 気象解析ライブラリを作る上で非常に便利な技などを 紹介していますので併せて参照ください. TIPS は本チュートリアルの最後尾にまとめて掲載しています.

なお, 紹介しているコードは上記モジュールの一部を チュートリアル用に簡略化したものです.

3 モジュールの仕様を決める

まずモジュールの仕様を決めます. ここでは以下の仕様のモジュールを作成することにします.

モジュールの空間名は NumRu::Derivative とします. NumRu は NUMericalRUby の略で, 数値計算や可視化関連のライブラリを 入れるモジュール空間です.

4 NArray 版の作成

4.1 メソッドを定義

仕様にしたがってメソッドを定義します. 以下は NumRu::Derivative の 0 版です. ユーザが呼び出すメソッドは cderiv です. 与えられた配列の両端の境界を 1 グリッド拡張して値を線形補完し(b_expand_linear_ext), ついでデータ配列, 軸配列の双方の中央差分を計算します(cdiff). 最後にそれらの商を返します.

derivative_ver0.rb:

TIPS

例外処理 - 引数のチェック / 配列の次元の指定方法 / NArray 同士の 2 項演算

4.2 テスト

上記で定義したメソッドをテストします. テストに際して以下の点をチェックすることにします.

以下にテストスクリプトを示します. 実際の製作過程では一項目ごとにテストを定義, 確認する 方が確実です. まとめて実行するとどこにバグが仕込まれているのかわからな くなってしまいます.

test_derivative_ver0.rb:

上記スクリプトをコマンドラインから実行するとテスト結果が出力されます. ただし, モジュールを定義したファイル(derivative_ver0.rb) がカレントディレクトリになくてはなりません.

% ruby test_derivative_ver0.rb 

**** equally spaced grid ****
**** single-D ****
dfdx - analytic (except boundary): [mean] 0.03922348996 [max] 0.06451071621
dfdx - analytic (except boundary): [mean] 0.0100170063  [max] 0.01636835692
error change from nx=11->21: 0.2553828409
**** multi-D ****
0.0
test exception successful

**** non-uniform grid ****
"x(11):"
NArray.float(11):
[ 6.28319, 5.68526, 5.14424, 4.6547, 4.21175, 3.81094, 3.44829, 3.12014, ... ]
"x(21):"
NArray.float(21):
[ 6.28319, 5.97675, 5.68526, 5.40799, 5.14424, 4.89335, 4.6547, 4.42769, ... ]
dfdx - analytic (except boundary): [mean] 0.01898250849 [max] 0.03216015526
dfdx - analytic (except boundary): [mean] 0.004934445822        [max] 0.01194461777
error change from nx=11->21: 0.2599469835

上記の結果から, 格子点数が倍になると誤差は約 0.25 倍となり, 定義した微分メソッドは 基本的に 2 次精度であることがわかります. 多次元配列に対する処理も成功しており, 引数のチェックも正しくなされています. なぜか不等間隔格子配列に対する結果の方が誤差は小さいですが, たまたま与えた配列が相性が良かったと考えられます.

以上の結果からテストは成功です. 新しいメソッドを付け加える場合はその都度テストを追加していきましょう. また既存のメソッドに手を加えた場合 あわせてテストも書き換えます. テストは頻繁に, かつ慎重に積み重ねることが重要です.

モジュールファイルにテストを埋め込む / ユニットテストフレームワーク を用いたテスト

4.3 ドキュメント記述

定義したメソッドのテストが成功したらドキュメントを書きましょう. Ruby ではスクリプトファイルの中に埋め込む事を意図して定義 されたドキュメントフォーマット RD が存在します. 今回もスクリプトに RD でドキュメントを埋め込んでしまいます. RD の書式についてはドキュメンテーションリンク集 の RD の項のリンクを参照ください.

以下にドキュメントを埋め込んだソースを示します. ここでは英語で書いてますが, 日本語でも構いません. (Ruby ユーザの間では英語のドキュメントを用意しておくというのが通例のようです.)

derivative_ver1.rb:

コツは一つのメソッドが完成するごとに記述していくことです. ちまちまと少しずつ付け足してます. メソッドの仕様を変更するたびに更新していくのも重要です. モジュール全体が完成してから書こうと思うとあとで大変な目に会います. (更新内容を忘れてしまうという以上にモチベーションが上がりません.)

以上をもって NumRu::Derivative ver 1 の完成です. 今後は新しいメソッドを追加したリ, 既存のメソッドを改良したりしてモジュールのブラッシュアップを図ります. ( 製品版の微分演算モジュールでは cderiv に加えて 3 点利用二次精度差分(threepoint_O2nd_deriv) が 定義されています. ) その都度, 上記の工程( メソッド定義/改良 --> テスト --> ドキュメント記述) を繰り返します.

5 GPhys 版 の作成

以上で, 微分メソッドはできました. しかし, このままでは, 座標とデータは 一まとまりになっておらず, 単位も名前もありません. そこで GPhys 版を作 りましょう. GPhys オブジェクトを対象に微分を行い, 結果もGPhysとして返 します. 微分すると一般に単位も変わりますので, 正しく自動更新される ようにしましょう. こうすることで, 自己記述的なデータから自己記述的な 計算結果を一行で得ることが出来るようになります. すると, あと 1-2 行で 図にすることも, ファイルに自己記述的な形で出力することもできますので, とても便利です.

実現方法としては, (1) GPhys にメソッドを追加する, (2) モジュール 関数にする, の2つが考えられますが, ここでは後者を採用します. 前者のほうが自然なのですが, 無闇に自分でメソッド追加をすると, GPhys 備え付けのメソッドと自前で追加したメソッドの区別が難しく なりがちという短所もあります. なお, これは GPhys 本家に取り込んで 貰いたいというのが出来ましたら, 積極的に作者まで連絡してください. (連絡先は電脳 Ruby のメーリングリストへ)

5.1 メソッド定義

作成した NumRu::Derivative を利用して GPhys 版 Derivative を作成します. 微分演算のエンジンは NumRu::Derivative にお任せにして, GPhys 版では物理量の単位と軸名などの属性に留意します.

gphys_derivative_ver0.rb:

微分演算部分は 22 行目だけです. 特に 27 行目以降はデータの単位と long_name 属性のケア を行っています.

TIPS

GPhys の作り方

5.2 テスト(解析解を用いて)

GPhys 版のテストでは以下の点に注意します. 計算誤差等のチェックは NArray 版の方で行っていますので ここではしません.

以下にテストスクリプトを示します. 長いですが, 解析階も含め自己記述的な GPhys データをゼロから作り出して るからです. 最初から自己記述的な NetCDF や GrADS ファイル中のデータ を読んで微分するだけなら, 劇的に短くなります(下記の 4.3.3 節 参照) .

test_gphys_derivative_ver0.rb:

上記スクリプトをコマンドラインから実行するとテスト結果が出力されます. derivative_ver0.rb, gphys_derivative_ver0.rb がカレントディレクトリになくてはなりません. 最初の行で警告が出ますが, gphys-0.3.5 に NumRu::Derivative および GPhys::Derivative が既に取り込まれているためです.

% ruby test_gphys_derivative_ver0.rb 

./derivative_ver0.rb:8: warning: already initialized constant LINEAR_EXT
******** dimname == 'lat' ********
dfdx - kaiseki_kai (except boundary): 0.01307030336     0.1701432467
**** check attribute ****
<attr-name>    <before>       <after>
name           t              dt_dlat
long_name      temperature    d_temperature_d_latitude
units          K              K.rad-1
******** dim == 0(lon) ********
dfdx - kaiseki_kai (except boundary): 0.01806355861     0.07927459478
**** check attribute ****
<attr-name>    <before>       <after>
name           t              dt_dlon
long_name      temperature    d_temperature_d_longitude
units          K              K.rad-1
******** dim == -1(z) ********
dfdx - kaiseki_kai (except boundary): 3.220621713e-08   7.748603821e-07
**** check attribute ****
<attr-name>    <before>       <after>
name           t              dt_dz
long_name      temperature    d_temperature_d_altitude
units          K              K.m-1

5.3 テスト(実際のデータを用いて)

ここでは GPhys パッケージに付随するテストデータ(T.jan.nc)を用いて テストをしましょう. このデータはNCEP 再解析による, 全球の気温の1月の気候値です. テストデータの詳細は GPhys のチュートリアル:本チュートリアルで用いるデータ をご覧ください.

以下テストスクリプトです.

test_gphys_derivative_ver1.rb:

01  require "numru/gphys"
02  require "gphys_derivative_ver0"
03  
04  include NumRu
05  
06  t = GPhys::IO.open('T.jan.nc','T')      
07  dtdx = GPhys::Derivative.cderiv(t,'lon')
08  newfile = NetCDF.create('dtdx.jan.nc')   
09  GPhys::IO.write(newfile, dtdx)          
10  newfile.close                           

T.jan.nc 中の変数 "T" の経度微分を dtdx.jan.nc というファイルに保存しています. 微分がたった 1 行(05 行目)であることに注目です. これだけで 温度場の経度微分が計算できます. 上記のスクリプトを実行して生成された netCDF ファイルの中身を見てみましょう.

% ruby test_gphys_derivative_ver1.rb
% ncdump -h dtdx.jan.nc |less

tcdf dtdx.jan {
dimensions:
            lon = 36 ;
            lat = 19 ;
            level = 9 ;
variables:
            float lon(lon) ;
                    lon:long_name = "Longitude" ;
                    lon:units = "degrees_east" ;
                    lon:actual_range = 0.f, 357.5f ;
            float lat(lat) ;
                    lat:long_name = "Latitude" ;
                    lat:units = "degrees_north" ;
                    lat:actual_range = 90.f, -90.f ;
            float level(level) ;
                    level:GRIB_id = 100s ;
                    level:positive = "down" ;
                    level:long_name = "Level" ;
                    level:units = "millibar" ;
                    level:actual_range = 10.f, 1000.f ;
                    level:GRIB_name = "hPa" ;
            double dT_dlon(level, lat, lon) ;
                    dT_dlon:parent_stat = "Other\n",
        "-" ;
                    dT_dlon:var_desc = "Air Temperature\n",
        "A" ;
                    dT_dlon:missing_value = -9.96921e+36f ;
                    dT_dlon:long_name = "d_Temperature_d_Longitude" ;
                    dT_dlon:least_significant_digit = 0s ;
                    dT_dlon:scale_factor = 1.f ;
                    dT_dlon:units = "degC.degrees_east-1" ;
                    dT_dlon:level_desc = "Multiple levels\n",
        "F" ;
                    dT_dlon:actual_range = -72.66724f, -24.35f ;
                    dT_dlon:dataset = "CDC Derived NCEP Reanalysis Products\n",
        "AC" ;
                    dT_dlon:add_offset = 0.f ;
                    dT_dlon:precision = 0s ;
                    dT_dlon:statistic = "Long Term Mean\n",
        "C" ;

// global attributes:
           :history = "2005-03-11 16:14:03 JST daktu32> NumRu::GPhys::NetCDF_IO.write dT_dlon" ;
}

更に gphys 付属の gpview でコンター図を見てみましょう.

% gpview dtdx.jan.nc@dT_dlon

このように, 単位などの属性のケアも 1 行でできるようになりました.

5.4 ドキュメント記述

NArray 版同様にドキュメントをつけて完成です.

gphys_derivative_ver1.rb:

6 完成

以上で微分演算モジュールが完成しました. 実際に GPhys::EP_Flux の中で用いられている微分演算モジュールは 上記の cderiv に加えて不等間隔格子に対して 2 次精度の差分を与える メソッド(threepoint_O2_deriv) が加わっていますが, 基本的には上記で示したとおりです.

より複雑なライブラリの作成も基本的にはここで紹介した流れで開発を進めていきます. 意外と簡単だと思われたのではないでしょうか? 自分でもできそうだ, と感じた人はぜひ新たなライブラリを作成してみてください. 繰り返しになりますが, GPhys 本家に取り込んでもらいたいものができましたら 積極的に作者まで連絡してください. また, 完成していなくてもこんなもの作ってます, という PR も大歓迎です.


7 (付録) TIPS

ここでは本文で紹介したライブラリコード中で登場する小技(TIPS)について 解説をします. これらの TIPS は気象解析ライブラリを作成知る上で非常に 便利ですので, ぜひとも参考にしてください.

7.1 例外処理 - 引数のチェック

4.2.1 節で提示したコードの各メソッド定義部の先頭では引数のチェックが行われています. これは想定外の引数が渡された場合は「例外」を上げるというものです. Ruby では例外処理が言語仕様に組み込まれています. 組み込みメソッド raise は例外を発生させる命令で, 以下の書式にしたがって記述します.

raise 例外クラス, 警告文(String)

ここで例外クラスには引数に関する例外クラス(ArgumentError)などを入れます. ちなみに例外クラスは省略可能で

raise 警告文(String)

のようにも書けます. この時は RuntimeError が上がります.

とある条件の時に例外を挙げる場合は

if 条件式

  raise HogeError, "herohero"

end

raise HogeError, "herohero" if 条件式

のように書きます. より詳しく知りたい方はこちら をご覧ください.

7.2 配列の次元の指定方法

4.2.1 節で提示したコードでは以下のような配列の次元の指定の仕方が多用されています. これは配列 zの dim 次元目の先頭の値を取り出すものです.

42       val0  = z[*([true]*dim +  [0] + [false])]      # 配列の先頭の値

これは以下の表記と同じです.

val0  = z[true, tru, true, ..., 0, false]      
          ^^^^^^^^^^^^^^^
           true が dim 個

ここで true はその次元の配列要素全てを表しています. また false は任意個の true を表す「rubber(ゴム)次元」です. Array どうしの足し算は両辺をくっつけ, 掛け算は個数分の繰り返しであること, そして最後の引数が * のついた配列である場合, その中身を展開することを利用してます. ( [dennou-ruby:002026] 参照)

ここで Ruby の配列要素の指定の仕方について触れておきます. Ruby の配列は先頭から 0, 1, 2... と数えます. 一方後ろから数える場合は -1, -2, .. となります.

ary = [2, 3, 5, 9, 10]
p ary[0]  #=> 2
p ary[-1] #=> 10  
p ary[1..-2] #=> [3, 5, 9]

7.3 NArray 同士の 2 項演算

NArray 同士の 2 項演算では長さ1の次元は任意の長さに自動的に拡張されます. 具体例をお見せしましょう. 以下は [3,3] の配列と [3,1] の配列の和の 演算です. この場合 [3,1] の配列の 2 次元目は自動的に長さ 3 に拡張 されます.

Narray[ [0,1,2],
        [3,4,5],     +   NArray[ [ 0,100,200 ] ]
        [6,7,8]]

---> # 内部で自動的に拡張

Narray[ [0,1,2],         NArray[ [ 0,100,200 ]
        [3,4,5],     +           [ 0,100,200 ]
        [6,7,8]]                 [ 0,100,200 ] ]

---> # 結果
Narray[ [0,101,202],     
        [3,104,205],     
        [6,107,208]]     

7.4 モジュールファイルにテストを埋め込む

4.2.2 節で示したテストスクリプトはモジュールを定義したファイルと独立にしましたが, 短いスクリプトならばモジュール本体のファイルに埋め込む方がよいでしょう. 無駄にファイルを増やさずに済みますし, モジュールの変更に応じてテストを書き 換えたり実行するのが簡単です.

テストを埋め込むにはテストコードを

if $0 == __FILE__

end

で挟みます. $0 は ruby インタプリタによって実行されている ファイルを表す変数, __FILE__ は自身のスクリプトファイルを表す定数です. こうするとモジュールファイル自身を実行した場合のみ, テストコードが実行されます. モジュールファイルを require してもテスト部分は無視されます.

7.5 ユニットテストフレームワーク を用いたテスト

Ruby にはユニットテストを支援するためのクラス RubyUnit(TEST::Unit) が存在します. 元々外部ライブラリだったのですが, Ruby 1.8 から標準添付ライブラリとして TEST::Unit というクラスになりました. 個人的な感想として数値計算モジュールのテストにはあまり向かない (ex. 精度の確認など) と思っているのですが, 興味のある方は利用して見てください. そして数値計算向けのテストの方法を編み出した方, 積極的に申し出てください.

7.6 VArray/GPhys の作り方

4.3.2 節で示したコードの 24-25 行目では微分結果の NArray を元に VArray, そして GPhys を構成しています.

24         v_dgpdx = VArray.new(n_dgpdx, v_data, name)            # 微分結果の VArray を構成
25         g_dgpdx = GPhys.new( gp.grid_copy, v_dgpdx )                        GPhys  を構成

VArray は基本的にデータ本体の NArray に, そのデータの名前と属性を加 えた配列です. VArray クラスのインスタンスを生成するには

VArray.new( データ本体(NArray), 属性(Attribute, Hash, or VArray; 省略可, 変数名(String; 省略可) )

のようにします. 上記のように第二引数に VArray を与えると, その属 性がコピーされて, 新たに生成される VArray オブジェクト に引き継がれます. 属性を陽に定義したいときは Hash を 使うと良いでしょう. なお, NumRu::Attributeはgphysのパッケージ内 で定義されています(解説は略). 第二, 第三引 数ともに省略可能で, 変数名は "noname" となります。nilを与えても 省略したことになります. 上記の例ではデータ本体を微分結果の NArray オブジェクト, 属性は微分前の属性を引き継ぎ新たに 生成した name を用いています.

GPhys はデータ配列(VArray) と格子情報(Grid) からなります. GPhys クラスのインスタンス生成するには

GPhys.new( 格子情報(Grid), データ配列( VArray or VArrayComposit ) )

のようにします. 先の例(25 行目)では 24 行目で構成した VArray をデータ配列して 元の GPhys オブジェクトのグリッドをコピーして( gp.grid_copy )新たな GPhys オブジェクトを生成しています.