はじめに

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

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

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

概要編
Dataset編
DataArray編 (本記事)

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

フォーマットの変換

既存の`xarray.DataArray`と同じ形の空の`xarray.DataArray`を作成する

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

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

緯度と経度を指定して空の`xarray.DataArray`を作成する

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`に変換する

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

変数・座標軸

座標軸

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

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

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

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

座標軸の名前を変更する

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

時間軸の扱い

ds.time # xarray.DataArray

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

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

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

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

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

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`の変数名自体を変更する

da.name = "SST"

配列操作

範囲選択、切り出し

インデックスによるスライス

numpy.ndarrayと同様

da[0,::-1,10:20] # xarray.DataArray

緯度経度を指定して、最も近い格子点の値を取得

直接緯度や経度を指定して取得することができる。

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

基本的にこの方法をとる必要がないが、インデックスの値を陽に求める古典的な方法は以下の通り。

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に変換する必要がある。

緯度経度など軸の範囲を指定してスライス

da.sel(
	time=slice(np.datetime64("2010-03-01"), np.datetime64("2011-02-28")),
	lat=slice(20, 50),
	lon=slice(120, 150))

特定の月のみを選択

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

座標軸の範囲を広げる

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で軸の名前を指定する。

da.mean(dim="time")

移動平均

dataarray.rollingを使う。pandasと同様。

da.rolling(time=3, center=True).mean() # 3時刻の移動平均。端にはnp.nanが入る。

時間間隔を粗くする（hourlyからdailyにする）

dataarray.resampleを使う。pandasと同様。

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

時間間隔を細かくする（線形内挿）

時間方向に内挿する例。monthlyデータからdailyにするときなど。

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に戻す（他にいいやり方があれば教えてください）

# 5km格子から1km格子に変換する(6x5)
da.values.repeat(6, axis=0).repeat(5, axis=1)

n*m格子ごとに平均して格子を粗くする

# 1km格子から5km格子に変換する(6x5)
da.coarsen(latitude=6, longitude=5).mean()

デフォルトではboundary="exact"なので、割り切れないとエラーになる。
da.coarsen()する前にda.sel()またはスライスして、開始・終了格子の位置を合わせておくとよい。

n*m格子ごとに、ある値を持つ格子が何格子あるか数える

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)

欠損値、マスク

マスクの取得

da.isnull() # NaNならTrue, それ以外はFalse
da.notnull() # NaNならFalse, それ以外はTrue

マスクする

da.where(mask) # maskがFalseの所の値がNaNになる

maskがxarray.DataArrayの場合、daと軸情報が同じである必要がある。一致しない場合は、mask.valuesによりnumpy.ndarrayに変換してから渡せばよい。

関数を適用してmapする

def myfunc(v):
    // do something
    return ret
    
xr.apply_ufunc(np.vectorize(myfunc), da)

描画 (クイックルック)

xarray.DataArray.plotを使うと、簡単に変数を描画できる（pandas.DataFrame.plotと同様）。

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

折れ線グラフ

da[:,0,0,0].plot()

2次元平面図(contour, shade)

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

2次元平面の矢印(quiver)

2変数必要になるので、xarray.Datasetにする必要がある

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")