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

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

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

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

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

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


## インストール
xgradsはcondaやconda-forgeでは配布されていないので、pipまたはgithub経由でインストールすることになります。Anacondaを使って仮想環境を管理したい場合は、同じライブラリをcondaとpipで重複して管理すると[面倒なことになる]("https://onoz000.hatenablog.com/entry/2018/02/11/142347")ので、必要な依存ライブラリをできるだけcondaでインストールしてから、xgradsだけをpipでインストールすることにします。

### 参考
- [CondaとPip (Python公式)]("https://www.python.jp/install/anaconda/pip_and_conda.html")
- [conda と pip でのパッケージ管理について]("https://qiita.com/ynakayama/items/29efebeb38604d10acef")

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


```zsh
conda create -n xgrads_test xarray dask numpy cartopy pyproj
```

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

```zsh
conda activate xgrads_test
conda install jupyter jupyterlab matplotlib
```

満を持してpip install。

```zsh
pip install xgrads
```

## GrADSファイルを見てみる


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

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

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

```python
from xgrads import CtlDescriptor, open_CtlDataset
import xarray as xr

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

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

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

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

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

```python
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"])`などを用いることで確認できます。

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

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

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