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

# はじめに
[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/items/586960f08c8b9fc30542)
2. [Dataset編](https://qiita.com/wm-ytakano/items/ddeecf2b43a4d1b200f2)
3. DataArray編 (本記事)

第3部となる本記事では、`xarray.DataArray`に焦点をあて、実際のデータ解析でよく使う機能を紹介する。

# フォーマットの変換
## 既存の`xarray.DataArray`と同じ形の空の`xarray.DataArray`を作成する

作成した`xarray.DataArray`の`values`に`numpy.ndarray`を代入することで、`numpy.ndarray`に軸情報を加えた`xarray.DataArray`を作成することができる。

```Python
xarray.full_like(other, fill_value, dtype=None)
xarray.zeros_like(other, dtype=None)
xarray.ones_like(other, dtype=None)
```

## 緯度と経度を指定して空の`xarray.DataArray`を作成する
```Python
latitudes = np.linspace(-90, 90, 181)
longitudes = np.linspace(-180, 180, 361)

empty_data = np.empty((len(latitudes), len(longitudes)))
empty_data[:] = np.nan

coords = {
    "latitude": (["latitude"], latitudes),
    "longitude": (["longitude"], longitudes),
}

empty_da = xr.DataArray(empty_data, dims=("latitude", "longitude"), coords=coords) 
```

## `xarray.Dataset`に変換する

```python
da.to_dataset(name="変数名")
```

# 変数・座標軸
## 座標軸

### `xarray.DataArray`から座標軸情報を取得する

```python
da.coords
da.coords["time"] # da.coords.timeではアクセスできない
```

### xarray.DataArrayの座標軸の情報を変更する（上書きする・書き換える）

```python
da.assign_coords({"time": time})
```

### 座標軸の名前を変更する

```python
da_new = da.rename({"latitude": "lat", "longitude": "lon"})
```

## 時間軸の扱い

```python
ds.time # xarray.DataArray
```

`xarray.DataArray`の形なら`.dt`アクセサを使って、`datetime.datetime`と同様の扱いができる

```Python
ds.time.dt.year # xarray.DataArray <int>
ds.time.dt.strftime("%Y/%m/%d") # xarray.DataArray <str>
```

`.values`にするとnp.datetime64型の配列になる

```Python
ds.time.values # numpy.ndarray(dtype='datetime64[ns]')
```

これを`datetime.datetime`型に変換するにはpandasの`pd.to_datetime`を使うのが早い

```Python
import pandas as pd

# 時刻が入ったxarray.DataArray
valid_time = ds["valid_time"] 
# numpy.datetime64型のある時刻を取り出す
npdatetime = valid_time[0].values
# pandas.Timestamp型に変換
pddatetime = pd.to_datetime(npdatetime)
# pandas.Timestampはdatetime.datetimeをPandasでラップしたものなので、datetime.datetimeと同じように使える
pddatetime.strftime("%Y%m%d") # 20220301
# 本当のdatetime.datetimeが欲しければ、以下のようにする
pddatetime.to_pydatetime() # datetime.datetime
```

## 変数

### `xarray.DataArray`の変数名自体を変更する

```python
da.name = "SST"
```

# 配列操作

## 範囲選択、切り出し
### インデックスによるスライス
`numpy.ndarray`と同様
```python
da[0,::-1,10:20] # xarray.DataArray
```

### 緯度経度を指定して、最も近い格子点の値を取得
直接緯度や経度を指定して取得することができる。

```python
da.sel(lon=90, lat=0, method='nearest')
```

基本的にこの方法をとる必要がないが、インデックスの値を陽に求める古典的な方法は以下の通り。
```python
import numpy as np
i = np.argmin(np.abs(ds.lon.values - 120))
j = np.argmin(np.abs(ds.lat.values - 20))
ds.sst[:, j, i]
```
`xarray.DataArray`(`ds.lon`, `ds.lat`)のままだと、`np.argmin`がエラーになるので、`ds.lon.values`として`numpy.ndarray`に変換する必要がある。

### 緯度経度など軸の範囲を指定してスライス
```python
da.sel(
	time=slice(np.datetime64("2010-03-01"), np.datetime64("2011-02-28")),
	lat=slice(20, 50),
	lon=slice(120, 150))
```

### 特定の月のみを選択

```python
da.sel(time=da['time.month']==1) # 1月のみを選択
```

### 座標軸の範囲を広げる

```python
da_new = da.reindex(
    {"latitude": new_latitudes, "longitude": new_longitudes}, 
    method="nearest", 
    fill_value=0)
```

`method="nearest"`を指定しないと、新しい座標が古い座標と少しでもずれていると値がfill_valueになってしまう。
`fill_value`で拡張された領域を埋める値を指定できる。デフォルトはNaN。

## 時間方向の解析

### 時間平均
`np.ndarray`の場合は`axis`で指定するが、`xarray.DataArray`の場合は`dim`で軸の名前を指定する。

```python
da.mean(dim="time")
```

### 移動平均
`dataarray.rolling`を使う。pandasと同様。
```python
da.rolling(time=3, center=True).mean() # 3時刻の移動平均。端にはnp.nanが入る。
```


### 時間間隔を粗くする（hourlyからdailyにする）
`dataarray.resample`を使う。pandasと同様。

```python
da.resample(time="1D").mean() # 平均する
da.resample(time="1D").count() # データの個数(=24)、正しく平均等の処理がされているかの確認に使える
```

### 時間間隔を細かくする（線形内挿）
時間方向に内挿する例。monthlyデータからdailyにするときなど。
```python
time = np.arange(
 	np.datetime64("2000-01-01"), 
 	np.datetime64("2001-01-01"), 
 	np.timedelta64(1, "D"))
da.interp({"time": time}, method="linear")
```

## 空間方向の解析

### n*m格子ごとに同じ値を繰り返して格子を細かくする

`numpy.ndarray.repeat`で値を作ったあと、軸を設定して`xarray.DataArray`に戻す（他にいいやり方があれば教えてください）
```python
# 5km格子から1km格子に変換する(6x5)
da.values.repeat(6, axis=0).repeat(5, axis=1)
```

### n*m格子ごとに平均して格子を粗くする
```Python
# 1km格子から5km格子に変換する(6x5)
da.coarsen(latitude=6, longitude=5).mean()
```
デフォルトでは`boundary="exact"`なので、割り切れないとエラーになる。
`da.coarsen()`する前に`da.sel()`またはスライスして、開始・終了格子の位置を合わせておくとよい。

### n*m格子ごとに、ある値を持つ格子が何格子あるか数える
```python
def count_mesh(array, axis=None, **kwargs):
    return np.sum(array == 1, axis=axis)
da_count_10x10 = da.coarsen(latitude=10, longitude=10).reduce(count_mesh)
```

## 欠損値、マスク
### マスクの取得
```python
da.isnull() # NaNならTrue, それ以外はFalse
da.notnull() # NaNならFalse, それ以外はTrue
```

### マスクする
```python
da.where(mask) # maskがFalseの所の値がNaNになる
```
`mask`が`xarray.DataArray`の場合、`da`と軸情報が同じである必要がある。一致しない場合は、`mask.values`により`numpy.ndarray`に変換してから渡せばよい。

## 関数を適用してmapする
```python
def myfunc(v):
    // do something
    return ret
    
xr.apply_ufunc(np.vectorize(myfunc), da)
```

# 描画 (クイックルック)
`xarray.DataArray.plot`を使うと、簡単に変数を描画できる（`pandas.DataFrame.plot`と同様）。

なお、以前書いた記事と内容が重複するところがある。

https://qiita.com/wm-ytakano/items/01a9b08885439d0adfb4

## 折れ線グラフ
```python
da[:,0,0,0].plot()
```

## 2次元平面図(contour, shade)

```python
# 変数はあらかじめ(緯度経度の)2次元にしておく必要がある
da[0,0,:,:].plot.contour() # 等値線
da[0,0,:,:].plot.contourf() # 色の陰影
da[0,0,:,:].plot.pcolormesh() # 色の陰影
```

## 2次元平面の矢印(quiver)
2変数必要になるので、`xarray.Dataset`にする必要がある

```python
rad = np.radians(da_dir8.astype(np.float32) * 45) # 8方位→ラジアン
ds_uv = xr.Dataset({"U": np.sin(rad), "V": np.cos(rad)})
ds_uv.plot.quiver(x="lon", y="lat", u="U", v="V")
```