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

はじめに

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

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

以前はGrADSを使ってクイックルックをすることが多かったのですが、Pythonに移行する流れになりつつあります。この記事では、xarray.DataArray.plotを使ったクイックルックの方法について説明します。

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

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つがインポートされているものとします。

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

Xarrayを用いたデータの読み込みとデータ型

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

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

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

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

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

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

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

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

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

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

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

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）は以下のように分類されます（日本語訳がないのでそのまま英語で書きます）。

• 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次元になるようにスライスしておく必要があります。

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

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

# 描画する
da.plot.pcolormesh() 

細かな注意点は以下の通りです。

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

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

da.plot.pcolormesh(center=False) 

さらに色の塗り方を調整したい場合は、cmapやvmin, vmaxを指定するとよいでしょう。その他オプション引数の詳細については公式ドキュメントをご参照下さい。

クイックルック（海岸線あり）

(論理)行数が2行増えますが、xarray.DataArray.plotで海岸線を描くことも可能です。

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='--')

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

p = da.plot(
    subplot_kws={"projection": ccrs.NorthPolarStereo(central_longitude=140)},
    transform=ccrs.PlateCarree())

謝辞

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