Xarrayによるモダンなデータ解析(2) Dataset編

はじめに

Xarrayは大気海洋分野のデータ解析に用いるPythonパッケージのデファクトスタンダートとなっている。しかし、日本語の情報があまりまとまっておらず、誰もが十分に使いこなせているわけではない。

以前「PythonでnetCDFを読み込む際のnetCDF4とxarrayの違い」という記事を書いたが、特にnetCDF4-pythonパッケージから移行した人は昔ながらの最低限の機能を使って書いており、最近の便利な機能を知らないことも多いと思われる。

そこで、以下の3部構成で"モダンな"データ解析方法について解説する。

1. 概要編
2. Dataset編 (本記事)
3. DataArray編

本記事では、xarray.Datasetに焦点を当て、ファイルの読み書きやデータフォーマットの変換について説明する。

ファイル入出力

netCDF

1つのnetCDFを読み込む

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

複数のnetCDFを読み込む

まとめて一つのDatasetとして読み込むことができる

ds = xr.open_mfdataset("*.nc")

• 内部でDaskが使われているので、インストールが必須

  • $ conda install -c conda-forge dask

• Daskが使われているため、openしてもメモリには読み込まれる、計算結果は遅延評価される。メモリ上に載せるにはds.compute()やds.load()を使う

• 読み込み時にHDF5の以下のようなエラーが出ることがある

  HDF5-DIAG: Error detected in HDF5 (1.12.2) thread 1:
      #000: H5A.c line 528 in H5Aopen_by_name(): can't open attribute
          major: Attribute
          minor: Can't open object
  

  この場合は、以下のようにする

  from dask.distributed import Client
  c = Client(n_workers=1, threads_per_worker=1)
  

  参考：https://stackoverflow.com/questions/72821108/hdf5-warnings-when-accessing-xarray-dataset

netCDFに出力する

ds.to_netcdf("example.nc")

xarray.Datasetを年月別に分割してnetCDF出力する

xarray.Dataset.groupBy()を使う。（ただしこの方法だと、ファイルサイズが大きくなりすぎる場合、出力時にメモリを大量に消費してしまう）

for year, ds_year in ds.groupby("time.year"):
    for month, ds_month in ds_year.groupby("time.month"):
        print(year, month)
        ds_month.to_netcdf(f"{year}{month:02d}.nc")

to_netcdfの代わりにsave_mfdatasetを使う方法もある（未調査）
https://docs.xarray.dev/en/stable/generated/xarray.save_mfdataset.html

netCDF4出力時に圧縮する・チャンクを指定する

ds.to_netcdf(
    filename,
    encoding={'dswrf': {
        "zlib": True,
        "complevel": 1,
        "chunksizes": (16,16,78),
}})

• 変数(xarray.DataArray)ごとに圧縮やチャンクサイズを指定する必要がある
• complevelは0で非圧縮、1以上の数字が大きくなると圧縮率が高い
  • 個人的な経験だと、complevelは1でよい。complevelを2以上にしても圧縮にかかる時間が長くなる割に、あまりサイズは小さくならない。

RuntimeError: NetCDF: HDF error

ディスクが溢れて書き込めない時に出ることがある
参考： https://github.com/pydata/xarray/issues/2132

ValueError: All sources must be dask array object

Metpyを使っているとき、metpy.quantify()したものがdatasetに含まれる場合に出る。metpy.dequantify()する

GRIB2

一般的なテーブルを使っているものは読める

ds = xr.open_dataset("example.grib", engine="cfgrib") # xarray.Dataset

• cfgribのインストールが必要
  $ conda install -c conda-forge cfgrib

• 気象庁のGRIBファイルなど独自のテーブルを使っている場合、読めないことがある。その場合はwgrib2等でnetCDFに変換して読み込む。

GrADS

プレーンバイナリ(.grd)とコントロールファイル(.ctl)が存在する場合は、xgradsパッケージを使うことでxarray.Datasetとして読み込める。

condaにはないのでpipでインストールする
$ pip install xgrads

from xgrads import open_CtlDataset
ds = open_CtlDataset('test.ctl')

GTOOL3

MIROCやCOCOなどのGCMで使われているファイル形式。xgtool3というパッケージで読み込める

$ pip install git+https://github.com/k1bb/xgtool3

import xgtool3 as xgt3
path = 'your_data_dir/ATM/Ts'
file = xgt3.Gtool3(path)
da = file.open()

Xarrayにこだわらなければ、gtool3パッケージの方がPythonでは普及している。解説は https://isotope.iis.u-tokyo.ac.jp/~kanon/src/gtool4python.html を参照。

GeoTIFF

rioxarrayを使う。xarray.open_rasterio()はdeprecated

$ conda install -c conda-forge rioxarray

import rioxarray
da = rioxarray.open_rasterio(file_path, masked=True)

masked=Trueをつけないと、_FillValueの部分がマスクされない

Zarr

Zarrに出力する

ds.to_zarr("hoge.zarr")

Zarrのチャンクはds.chunk()で指定できる

ds.chunk({"longitude": 288, "latitude": 145, "time": 1}).to_zarr("hoge.zarr")

Xarrayで作成したZarrを読み込む

xr.open_zarr("hoge.zarr")

ファイル形式の変換

xarray.DataArray→xarray.Dataset

ds = xr.Dataset({
    "SST": da_sst,
    "Precip": da_precip,
})

複数のDataArrayをまとめてDatasetを作成するときに、時間軸がずれていると不都合がある。
その場合はda.assign_coords()で軸情報をあらかじめ書き換えて揃えておく。

numpy.ndarray→xarray.Dataset

ds = xr.Dataset(
   data_vars={
       "low": (["lat", "lon"], low, {"units": "percent", "long_name": "Low Cloud Cover"}),
       "middle": (["lat", "lon"], middle, {"units": "percent", "long_name": "Middle Cloud Cover"}),
       "high": (["lat", "lon"], high, {"units": "percent", "long_name": "High Cloud Cover"}),
   },
   coords={
       "lon": (["lon"], lon, {"units": "degrees_east", "long_name": "longitude"}),
       "lat": (["lat"], lat, {"units": "degrees_north", "long_name": "latitude"}),
   },
)

他のツールで使用することを考えると、netcdfのCF規約にしたがってunits属性とlong_name属性を加えておくとよい。

pandas.DataFrame→xarray.Dataset

時間・緯度・高度などの軸は、DataFrameの(マルチ)インデックスとして指定しておく必要がある

ds = xr.Dataset.from_dataframe(df)

xarray.Dataset→pandas.DataFrame

これはxarray.DataArrayにも適用可

ds.to_dataframe()

マルチインデックスの順番を指定できる

ds.to_dataframe(dim_order=["time","lat","lon"])

演算

Dataset同士を結合する

ds = xr.merge([ds_1, ds_2])

Dataset同士の四則演算

含まれるDataArrayの変数名が共通していて、共通する座標軸を持っている場合、まとめて計算できる。変数方向のループを書かなくてよいため、使えるとかなり便利。例えばDataset内の変数全てに対して、時間平均と標準偏差を計算し、それを用いて規格化するといったことが可能。

ds # (time, lat, lon)の形の変数が複数含まれているとする
ds_mean = ds.mean(dim="time") # (lat, lon)
ds_std = ds.std(dim="time") # (lat, lon)
ds_norm = (ds - ds_mean) / ds_std # (time, lat, lon)