-- GPhysレシピ --

格子/座標の取扱い

堀之内武 最終更新2012/12/10

GPhys を使ってちょっと凝ったことをすると,格子や座標を陽に扱うことがしばしば必要になるでしょう. この文書では GPhys で格子や座標をどう扱えるか,例で示していきます.

目次

準備

NCEP 再解析の等圧面高度 hgt.2010JanFeb.nc を適当な作業ディレクトリにダウンロードし,そこに cd します.

はじめに

まず,基礎となる GPhys の構造について説明します. GPhysオブジェクトは下図のように主変数(配列データ)と格子からなります.

gphys_structure_th.png

ここで,"has x" は x 個持つことを表します.has 1.. は1以上という意味です.

格子の中には,主変数の各次元に対応する座標軸があります. 座標軸は重複のない名前をもち,名前だけで次元を特定できます(名前が重複するようには作れないようになってます). 各座標軸は,それに沿った格子点の位置を保持する一次元配列を持ちます. 座標にがなじまない次元にも何か名前をつけ,位置もダミーでつけます (0,1,2,3,.. のインデックスなど).各座標軸は任意個補助変数を持てます.

例えば NetCDF や GRIB ファイルの中のある変数を GPhys::IO.open で「開く」際は, それぞれのファイル形式の自己記述性を活かして,座標までも芋づる式にたどって GPhys オブジェクトを構成します.NetCDF なら通常はすべての情報はファイルからとられるので 座標値もファイル中の変数となります(具体的にはVArrayNetCDFオブジェクトに格納される) *1. GRIB なら例えばヘッダー中の「開始位置」と「増分」より等間隔の1次元配列を生成したうえで保持する(VArrayとして)などとなります.

[*この段落は進んだ内容なので読み飛ばして構わない]: 格子は,座標軸に加え,任意個の「補助座標」(AssocCoords;一般に多次元)を持つこともできます.これは例えばこんなときに使います.「主変数が極座標格子点上の時系列という3次元配列で,座標軸は方位角θ,中心からの距離r,時刻 t である.r,θに加え,各格子点の直交直線座標 x, y の値も格納したいので,2次元の補助座標として x, y を定義する.」 補助座標は,それ自体が GPhys オブジェクトです.ただし,主変数と座標軸を共有するという縛りがあります.上の例でいえば,主変数が T(r,θ,t) と3次元,補助変数が x(r,θ), y(r,θ) と2次元で,x, y は同一であるといったことです.

座標値を取り替える

以上の各要素は,同じ GPhys オブジェクト内で「取り換え可能」です. ここでは,ファイルから開いた GPhys オブジェクトについて, ファイルの中身を書き換えないで,プログラム上で座標軸を取り替えて使う場合を紹介します.

次の例では,単位を hPa から Pa に変換し,名前を Level から pressure に付け替えたもので,鉛直座標を置き換えます. 比較のためもう一度同じのをオープンして並べて図示します. これ自体は特に意味のない例ですが,このあと例えば鉛直座標の単位が "Pa" であることを前提にした計算をしたいといった場合に役立ちます.

結果: replace_ax_th.png(クリックでフルサイズ表示)

コード解説: GPhys#axis メソッドは,引数(整数または文字列)で指定した座標 (Axis オブジェクト)を返します.Axis#pos メソッドは,その 座標の座標値を格納する1次元データ(VArray オブジェクト)を返しま す.よって,p = gp.axis(2).pos で得られる p は1次元の VArray オブジェクトです(今の場合,VArrayのサブクラスの VArrayNetCDF オブジェクトで,データ実体は NetCDF ファイル中にあります). これを単位変換し,長い名前(long_name)を付け替えて, gp.axis(2).set_pos(ppa) により,もとの (Axis オブジェクトにおいて座標値として与えます.set_pos は Axis オブジェクト内で置き換えを行うだけなので,もとのファイルには影響は及びません.

ちなみに,そもそも gp = GPhys::IO.open("air.2010JanFeb.nc", "air") では読み込み専用でファイルを開きますので,間違えて書こうとしても書けません. もしもファイルを書き込み可能な形で開きたければ,次のようにする必要があります:

file = NetCDF.open("air.2010JanFeb.nc", "a")   # "a"モードで書換可
gp  = GPhys::IO.open(file, "air")

補助座標(多次元)の導入

補助座標は座標変換などで威力を発揮します. GPhysオブジェクトに補助座標を導入するには GPhys#set_assoc_coords メソッドを使います. その例は 補間,座標変換 座標変換の節にあるので,参照してください.

座標や格子(やGPhysオブジェクトさえも)を一から構築

座標や格子も含め GPhys オブジェクトは一から(入力ファイル等を使わず)生成することができます.次はその例です.

結果: create_gphys_th.png(クリックでフルサイズ表示)

このように作った部品で,ファイルを開いて作ったGPhysの一部を置き換えることもできます (Axis#set_posレベルでも,ある座標軸全体でも,格子全体でも).


*1ただし,もしも座標変数がファイル中になければダミーが生成されます


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