GrADSファイルをXarrayに読み込む #Python

GrADS（およびGrADS形式のファイル）は大気・海洋分野では広く使われている描画ツールで、簡単に対話的な可視化を行うことができたり、さまざまな解析ができたりするのですが、普段それらの分野に携わっていないけどそういうデータを可視化したり解析したりしてみたいと思った時には、分野外の人間にとってはそれなりに尻込みしてしまうようなソフトウェア（およびファイル形式）でもあります。

できればPythonで扱いたいと思ってインターネットを探していたところ、以下のような選択肢を発見しました。

とにかく自分でいちからファイルを読みこむコードを書く
PyGradsやxgradsなどのライブラリを使う

前者はNumpyを使えば読み込めるといえば読み込めるのですが、データのさまざまな仕様を毎回きちんと理解して読み込むのに自信がなかったこと、手元にコントロールファイル (.ctl) があるならせっかくだしそれを使いたいなどの理由から、後者を選びました。PyGRADSはGrADSのコマンドをきちんと理解している人むけっぽい感じがしたのと、多次元データをいきなりXarrayに読み込んで使えるという点が良かったのでxgradsを使うことにしました。

私自身GrADSやそれをよく扱う分野に不慣れなので、つまりこの記事はそうした自分自身のためのインストールおよびチュートリアルの忘備録ということになります。

インストール

xgradsはcondaやconda-forgeでは配布されていないので、pipまたはgithub経由でインストールすることになります。Anacondaを使って仮想環境を管理したい場合は、同じライブラリをcondaとpipで重複して管理すると面倒なことになるので、必要な依存ライブラリをできるだけcondaでインストールしてから、xgradsだけをpipでインストールすることにします。

なので、最初にRequirementsをインストールします。ここではxgrads_testという名前の環境にしました。

conda create -n xgrads_test xarray dask numpy cartopy pyproj

のちのちデータを可視化しながらいじっていきたくなる気がするので、ついでにJupyterとかMatplotlibも。

conda activate xgrads_test
conda install jupyter jupyterlab matplotlib

満を持してpip install。

pip install xgrads

GrADSファイルを見てみる

インストールが終わったので試しになにかデータを見てみましょう。GrADSの開発元でもあるジョージ・メイソン大学の海洋陸面大気研究センター（COLA）の公式チュートリアルにサンプルデータ(example.tar.gz)があるのでこれを用います。example.tar.gzをダウンロード・展開するとexampleというフォルダができるのでその中で作業をします。

cd <exampleのあるディレクトリ>

CtlDexcriptor()でコントロールファイルを確認できます。

from xgrads import CtlDescriptor, open_CtlDataset
import xarray as xr

ctl = CtlDescriptor(file="model.ctl")
print(ctl)

GrADSファイルを読み込むにはopen_CtlDataset()を使います。

dset = open_CtlDataset("model.ctl")
print(dset)

コントロールファイルをパースしてGrADSファイルを読み込んでくれます。print(dset)の出力を見ればわかる通り、dsetの中身はxarray.Datasetです。

コントロールファイルが複数ある場合は、open_mfDataset()にワイルドカードを使って表したパスを渡してやることで読み込むことができます。

dset = open_mfDataset("*.ctl")

しばしば、1日ごとや1ヶ月ごとなどにデータが分かれており、（yyyymmdd.datのように）テンプレートに従ったファイル名が付けられていて、コントロールファイルの記述をパースすることでそれらを一括で読み込める仕様になっていることがあります（参考: http://cola.gmu.edu/grads/gadoc/templates.html）。

xgradsは大抵の場合をカバーしているのですが、たまにそのままでは読み込めないファイル名などがある（mm.dd.yyyyなど）ので、その場合はターミナルで一括リネームなどをしてから再度試してみるとうまく読み込めることがあります。

描画してみる

せっかくなので読み込んだデータを見てみましょう。COLAの公式チュートリアルでは最初にps（気圧）を描画しているので、ここでもそれをプロットしてみます（以下ではJupyterでの作業を想定しています）。
念の為ですが、psは地表面での気圧ですから、lev軸は存在しないことに注意が必要です。print(dset["ps"])などを用いることで確認できます。

fig = plt.figure()
ax = fig.add_subplot(projection = ccrs.Mercator(central_longitude = 180.0))

ax.coastlines()

dset["ps"][0,:,:].plot.contourf(
    ax = ax,
    transform = ccrs.PlateCarree(),
)

COLAの公式チュートリアルにある図と似た等圧線図が表示されると思います。
あるいは、海岸線が必要ないのであれば以下の一行で表示することができます。

dset["ps"][0,:,:].plot.contourf()

xgradsを使うことでGrADSファイルをいきなりxarrayに読み込むことができるため、xarray.DataArray.plotを使った多次元データのクイックルックと組み合わせることで、PythonでもGrADSデータの対話的な描画や、その先の解析がよりスムーズに行えるようになります。