はじめに
Xarrayは大気海洋分野のデータ解析に用いるPythonパッケージのデファクトスタンダートとなっている。しかし、日本語の情報があまりまとまっておらず、誰もが十分に使いこなせているわけではない。
以前「PythonでnetCDFを読み込む際のnetCDF4とxarrayの違い」という記事を書いたが、特にnetCDF4-pythonパッケージから移行した人は昔ながらの最低限の機能を使って書いており、最近の便利な機能を知らないことも多いと思われる。
そこで、以下の3部構成で"モダンな"データ解析方法について解説する。
第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")