Pythonのxarray.DataArray.plotがnetCDFのクイックルックに便利

# はじめに
大気海洋分野ではnetCDF形式のデータをよく用います。このデータは時空間方向の多次元配列を表現するバイナリデータであり、人間が把握するためには水平2次元の地図上に描画（俗に言う"お絵かき"）する必要があります。

研究などでデータ解析をしている最中は、手の込んだ図を描くよりもむしろ、さっとデータの中身を確認すること（クイックルック）が求められます。（どこかの先生が「最近の若い人は綺麗な図を描くようになったが、図を描くことばかりにこだわって肝心の解析をしなくなった」というような愚痴を言っていました。本当かは知りませんが。）

以前は[GrADS](http://cola.gmu.edu/grads/grads.php)を使ってクイックルックをすることが多かったのですが、Pythonに移行する流れになりつつあります。この記事では、`xarray.DataArray.plot`を使ったクイックルックの方法について説明します。

conda (miniforge)で環境構築しており、主要なパッケージのバージョンは以下の通りです。

```text
cartopy                   0.23.0
matplotlib                3.9.1
numpy                     2.0.1
python                    3.12.4
xarray                    2024.6.0
```

実行は、VSCode上で実行したJupyter Notebookを用いています。
以下の説明では`plt`,`xr`,`ccrs`の3つがインポートされているものとします。

```Python
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.crs as ccrs
```

# Xarrayを用いたデータの読み込みとデータ型
まず、データの読み込み方法と、Xarrayのデータ型についておさらいします。

例として"example.nc"という名前のMSM地上面の初期時刻のnetCDFを読み出すには次のようにします。読み出したファイルは、`xarray.Dataset`型として開かれます。

```Python
ds = xr.open_dataset("example.nc") # xarray.Dataset
ds
```

Jupyterではセルの最後の行に書かれた変数の内容が表示されるため、2行名に`ds`とだけ記述しています。するとデータセット`ds`の内容を確認できます。右側の📄や🗄️をクリックすることで、変数の属性や値を確認することができます。

![スクリーンショット 2024-09-12 17.03.47.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/829116/b7efa1fe-ce23-ff3e-a532-ab9d8f4cc9f3.png)


netCDFの中に含まれる"TMP_1D5maboveground"という多次元変数（1.5m気温）を読み出すことを考えます。`xarray.DataArray`型は軸の情報や属性(attribute)を含んだ型です。後で出てくるクイックルックで使います。

```Python
da = ds["TMP_1D5maboveground"] # xarray.DataArray
```

`numpy.ndarray`型にするには次のようにします。当然ながら軸情報も属性も含んでいないただの配列です。

```Python
# 以下どちらでも同じ
na = da.values # numpy.ndarray
```

`xarray.DataArray`, `numpy.ndarray`型はいずれもインデックスのスライスができます。スライスしても型は変わりません。以下の例では0番目の時刻を表す軸について、インデックスが0である最初の時刻を取り出しています。

```Python
da[0, :, :] # xarray.DataArray
na[0, :, :] # numpy.ndarray
```

`numpy.ndarray`になく`xarray.DataArray`にある特徴として、軸名を用いてスライスすることが可能です。`da.isel()`はインデックス、`da.sel()`は値で指定します。

```Python
da.isel(time=0)
da.sel(time="2024-04-01T00:00:00")
```

`da[0,:,:]`と`da.isel(time=0)`では結果は同じですが、後者の方がより可読性が高い書き方だと言えるでしょう。

# 2次元変数を描画する

## Matplotlib Application Interfaces

Matplotlibを用いた描画の方法・スタイル（Matplotlib Application Interfaces）は以下のように分類されます（日本語訳がないのでそのまま英語で書きます）。

https://matplotlib.org/stable/users/explain/figure/api_interfaces.html

- Native Matplotlib interfaces
    - explicit "Axes" interface (e.g., `ax.plot()`)
    - implicit "pyplot" interface：(e.g., `plt.plot()`)
- Third-party library "Data-object" interfaces： (e.g., `da.plot()`)

Matplotlib本体では2つの方法があり、一つは`fig`や`ax`を作成して描画するexplicit "Axes" interface、もう一つは`fig`や`ax`を用いずに`plt`のみで描画を行うimplicit "pyplot" interfaceです。explicit "Axes" interfaceの方が細かい図の設定が可能になる一方、記述量が多くなります。

一方、Third-party library "Data-object" interfacesは、XarrayやPandasなどの別のライブラリのデータオブジェクトに生えた描画メソッド`da.plot()`を用いる方法です。これは記述量が少なくクイックルックに有用です。

## クイックルック（海岸線なし）
Data-object interfacesである`xarray.DataArray.plot`を使うと、`import`文を除けば、ファイルオープン・変数の計算・描画の3行で図の描画が可能です。注意点として、`xarray.DataArray`はあらかじめ(緯度経度の)2次元になるようにスライスしておく必要があります。

```Python
# Datasetを開く
ds = xr.open_dataset("example.nc")

# DataArray（単位をKから℃に変換）
da = ds["TMP_1D5maboveground"].isel(time=0) - 273.15

# 描画する
da.plot.pcolormesh() 
```
![output.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/829116/79f698e2-cf92-5c1f-fe60-417a81b69aad.png)



細かな注意点は以下の通りです。
- JupyterではなくPythonスクリプトの場合`plt.show()`や`plt.savefig()`が必要です
- 何枚も図を描画する場合はメモリを解放するために`plt.close()`を実行してください
- ここでは2次元配列を`pcolormesh`で描画する例を示していますが、`contour`や`contourf`でも同様に描画できます

負の値が含まれている場合は、0を中心とした対称なカラーバーがデフォルトで設定されるようです。これを回避するには`center=False`を設定します。

```Python
da.plot.pcolormesh(center=False) 
```

![output2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/829116/8d5a6dd6-3f69-95b3-d347-c13197f81490.png)


さらに色の塗り方を調整したい場合は、`cmap`や`vmin`, `vmax`を指定するとよいでしょう。その他オプション引数の詳細については[公式ドキュメント](https://xarray.pydata.org/en/stable/generated/xarray.plot.pcolormesh.html)をご参照下さい。


## クイックルック（海岸線あり）
(論理)行数が2行増えますが、`xarray.DataArray.plot`で海岸線を描くことも可能です。

```Python
ds = xr.open_dataset("example.nc")
da = ds["TMP_1D5maboveground"].isel(time=0) - 273.15

# 座標系と投影法を指定して描画
p = da.plot(
    subplot_kws={"projection": ccrs.PlateCarree()},
    transform=ccrs.PlateCarree())

# 海岸線の描画
p.axes.coastlines()

# 緯線・経線の描画
p.axes.gridlines(
    crs=ccrs.PlateCarree(), 
    draw_labels=True, # 緯度経度ラベルを表示する
    linewidth=1, color='gray', alpha=0.5, linestyle='--')
```
![output4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/829116/771e0308-85f9-bd87-a063-31a94ab31206.png)


ここで`subplot_kws`の`projection`に投影法(PlateCarree, 等緯度経度)、transformにデータの座標系（こちらもPlateCarree）を指定しています。
`projection`を変えることで、別の投影法を選択することができます。例えば日本域の天気図でよく用いられるポーラーステレオ投影の場合は以下のように変更します。
```Python
p = da.plot(
    subplot_kws={"projection": ccrs.NorthPolarStereo(central_longitude=140)},
    transform=ccrs.PlateCarree())
```

![output5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/829116/bb809289-218e-ed2b-1c33-669bd38dc6fb.png)

# 謝辞
`center=False`のオプションについて、東京大学の佐藤瞭さんに情報提供頂き、記事に追記しました。ありがとうございました。