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

# はじめに
[Xarray](https://docs.xarray.dev/en/stable/index.html)は大気海洋分野のデータ解析に用いるPythonパッケージのデファクトスタンダートとなっている。しかし、日本語の情報があまりまとまっておらず、誰もが十分に使いこなせているわけではない。

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

そこで、以下の3部構成で"モダンな"データ解析方法について解説する。
1. [概要編](https://qiita.com/wm-ytakano/items/586960f08c8b9fc30542)
2. Dataset編 (本記事)
3. [DataArray編](https://qiita.com/wm-ytakano/items/564476a1e378a3bd7060)

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

# ファイル入出力

## netCDF

### 1つのnetCDFを読み込む
```python
ds = xr.open_dataset("example.nc") # xarray.Dataset
```

### 複数のnetCDFを読み込む
まとめて一つのDatasetとして読み込むことができる
```python
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
    ```
    この場合は、以下のようにする
    ```python
    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に出力する
```python
ds.to_netcdf("example.nc")
```

### xarray.Datasetを年月別に分割してnetCDF出力する
`xarray.Dataset.groupBy()`を使う。（ただしこの方法だと、ファイルサイズが大きくなりすぎる場合、出力時にメモリを大量に消費してしまう）
```python
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出力時に圧縮する・チャンクを指定する

https://qiita.com/wm-ytakano/items/73124ec9a3be8bc7b58c

### 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](https://github.com/miniufo/xgrads)パッケージを使うことでxarray.Datasetとして読み込める。

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

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

## GTOOL3
MIROCやCOCOなどのGCMで使われているファイル形式。[xgtool3](https://github.com/k1bb/xgtool3/tree/main)というパッケージで読み込める

`$ pip install git+https://github.com/k1bb/xgtool3`
```python
import xgtool3 as xgt3
path = 'your_data_dir/ATM/Ts'
file = xgt3.Gtool3(path)
da = file.open()
```

Xarrayにこだわらなければ、[gtool3](https://github.com/hyungjun/gtool3/)パッケージの方がPythonでは普及している。解説は以下の記事を参照。

https://qiita.com/kanonundgigue/items/82867eb43744bff53a1f 

## GeoTIFF
[rioxarray](https://corteva.github.io/rioxarray/stable/)を使う。`xarray.open_rasterio()`はdeprecated

`$ conda install -c conda-forge rioxarray`

```Python
import rioxarray
da = rioxarray.open_rasterio(file_path, masked=True)
```
`masked=True`をつけないと、_FillValueの部分がマスクされない

## Zarr
Zarrに出力する
```python
ds.to_zarr("hoge.zarr")
```

Zarrのチャンクは`ds.chunk()`で指定できる
```Python
ds.chunk({"longitude": 288, "latitude": 145, "time": 1}).to_zarr("hoge.zarr")
```

Xarrayで作成したZarrを読み込む
```Python
xr.open_zarr("hoge.zarr")
```
# ファイル形式の変換
## `xarray.DataArray`→`xarray.Dataset`
```python
ds = xr.Dataset({
    "SST": da_sst,
    "Precip": da_precip,
})
```
複数のDataArrayをまとめてDatasetを作成するときに、時間軸がずれていると不都合がある。
その場合は`da.assign_coords()`で軸情報をあらかじめ書き換えて揃えておく。

## `numpy.ndarray`→`xarray.Dataset`
```python
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の(マルチ)インデックスとして指定しておく必要がある
```python
ds = xr.Dataset.from_dataframe(df)
```

## `xarray.Dataset`→`pandas.DataFrame`
これは`xarray.DataArray`にも適用可
```python
ds.to_dataframe()
```

マルチインデックスの順番を指定できる
```python
ds.to_dataframe(dim_order=["time","lat","lon"])
```

# 演算
## Dataset同士を結合する
```python
ds = xr.merge([ds_1, ds_2])
```

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

```python
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)
```