xarrayでgradsの機能を代替

xarrayの詳しい説明はここ。

gradsでやっていたことをxarrayでできるようにする。gradsよりも軸の取り扱いが面倒なので、データを読み込んだ時点で必要に応じて軸データの変換を行う。

IO

gradsデータの読み込み

xgradsを使えば直接gradsデータをxarrayとして読み込める。ただし、

• 複雑な template
• yrev
• zrev
• sequential binary
  は苦手なので、その場合は１度netcdfデータに変換する。

from xgrads import open_CtlDataset

ds = open_CtlDataset("data.ctl")

netcdfの読み込み

ds=xr.open_dataset('input.nc')

DataArray, DataSetから変数を取り出す

以下の時間操作や空間操作では変数名の指定が抜けているが、基本的に変数を指定してデータを処理する。変数名を指定しない場合、データセット全てに対して同じ処理がかかる。

# safe method
u=ds['u']
# simple but not safe
u=ds.u
# multiple values
ds2=ds['u','v']

軸データの修正

# read data without decoding time
ds1=xr.open_dataset('temp_ens.nc',decode_times=False)
# rename time unit
ds1=ds1.rename({'time_counter':'time'})
# reset taxis by pandas
ds1["time"] = pd.date_range(start="2001-01-01",periods=276,freq="MS")

時間操作

時間平均

# monthly mean
da.groupby("time.month").mean("time")
# seasonal mean
da.groupby("time.season").mean("time")
# annual mean
annual = ds.t.groupby("time.year").mean()

時間平均（特定範囲）

# during specific periods (date)
tmean = ds.t.sel(time=slice("2001-01-01","2010-12-31")).mean("time")
# during specific periods (year only)
tmean = ds.t.sel(time=slice("2001","2010")).mean("time")
# during seasons
jja = ds.t.sel(time=ds.time.dt.month.isin([6,7,8])).mean("time")

アノマリー計算

# monthly climatology
clim = da.groupby("time.month").mean("time")
# monthly aomaly
anom = da.groupby("time.month") - clim

移動平均

# 3-month running mean, time centering, minimum 1 month data
da.rolling(time=3,center=True,min_periods=1).mean()

時間解像度変更

ds.resample(time='1D').sum() #daily sum
ds.resample(time='AS').mean() #annual mean

時間シフト相関

# da2 is 1 month before
corr = xr.corr(da1, da2.shift(time=1), dim="time")

空間操作

領域平均

import numpy as np

# latitudinal weight
weights = np.cos(np.deg2rad(ds.lat))

# average after weighting
mean_val = (
    ds.var.sel(lat=slice(10,20), lon=slice(120,150))
    .weighted(weights)
    .mean(("lat","lon"))
)

緯度経度サブセット

gradsだと

set lon 120 150
set lat -10 10

xarrayだと

da = ds.t.sel(
    lon=slice(120,150),
    lat=slice(-10,10)
)

zonal mean

#zonal mean
zonal = da.mean("lon")
#or
zonal = da.mean(("lon","time"))
#plot in latitude-level section
zonal.plot(
    x="lat",
    y="level",
    yincrease=False
)

zonal & meridional mean

u_eq = (
    ds.u
    .mean("lon")
    .sel(lat=slice(-5,5))
    .mean("lat")
)

厳密には緯度方向に重み付けをする。

import numpy as np

weights = np.cos(np.deg2rad(ds.lat))

u_eq = (
    ds.u.mean("lon")
    .weighted(weights)
    .mean("lat")
)

補間処理

以下は単純な緯度経度補間で球面上補間ではないので、xesmfを使ったほうがいい

ds.interp(lat=new_lat, lon=new_lon, method='linear')

log(p)で鉛直補間

da_log = da.assign_coords(lev=np.log(da.lev))

new_lev = np.log([925, 775, 600, 400])

da_interp = da_log.interp(lev=new_lev)
da_interp = da_interp.assign_coords(lev=np.exp(da_interp.lev))

リグリッド

xesmfを合わせて利用する。

import xarray as xr
import xesmf as xe
import numpy as np

ds_out = xr.Dataset(
    {
        "lat": (["lat"], np.arange(-89, 90, 2)),
        "lon": (["lon"], np.arange(0, 360, 2)),
    }
)

regridder = xe.Regridder(ds_in, ds_out, "bilinear")
ds_out = regridder(ds_in)

アンサンブル処理

アンサンブルデータ例

import xarray as xr
import numpy as np

da = xr.DataArray(
    np.random.rand(10,12,36,72),
    dims=("member","time","lat","lon")
)

ens_mean = da.mean("member")

ファイルが複数ある場合

以下で複数ファイルをアンサンブルメンバーのデータとして取り扱う

ds = xr.open_mfdataset(
    "run*.nc",
    concat_dim="member",
    combine="nested"
)

ens_mean = ds.var.mean("member")

統計値

# ensemble mean
ens_mean = ds.var.mean("member")
# variance
ens_var = ds.var.var("member")
# standard deviation
spread = ds.var.std("member")

条件抽出

ラベルによる抽出

ds.sel(time='2024-01-01', lat=35.0, lon=135.0)

マスク処理

# mask out nino34<1
da.sel(time=da["nino34"] > 1)

# mask out temp<308.15
ds.where(ds.temp > 308.15)

# mask out land=0
t_mask = ds.t.where(land==1)

描写

基本

da.plot()
ds.temp.isel(time=0).plot()

最も簡単なクイックルックコード

import xarray as xr
import matplotlib.pyplot as plt

# read data
ds = xr.open_dataset('sample_data.nc')
# display temp
ds.temp.isel(time=0).plot()
# display window
plt.show()

関数関係

回帰係数の算出(tregr)

# x: index, y: OLR
beta = xr.cov(x, y, dim="time") / xr.var(x, dim="time")

相関係数(tcorr)

xr.corr(x,y,dim="time")

鉛直積算(vint)

# precipitable water
g = 9.81
dp = q.level.diff("level")
q_mid = q.isel(level=slice(0,-1))
pw = (q_mid * dp).sum("level") / g

鉛直積算をかける前に、鉛直座標を必要に応じて反転させるには以下を使う。圧力に関して昇順にデータが並び替わる。

ds2=ds.sortby('plev',ascending=True)

回転, 発散(hcurl,hdivg)

残念ながらxarrayでなくてmetpyを使うのが簡単で確実。次元が合ってないと計算できないので注意。

import metpy.calc as mpcalc

# hcurl
curl = mpcalc.vorticity(u, v)
# hdivg
div = mpcalc.divergence(u, v)

metpyはunitがないとエラーになるので、xarrayでは以下でmetpyが読めるunitを与える

ds["u"].attrs["units"] = "m/s"
u = ds.u.metpy.quantify()

データセット

DataArrayは変数が一つのデータ、DataSetは変数が複数のデータ。

to_netcdfで簡単にnetcdfに出力できる。軸の情報が間違っていたり、違っているとxarrayで再度読み込んだときに計算に使えないので注意が必要。

da.to_netcdf('output.nc')

DataArray

da = xr.DataArray(
    np.random.rand(10,5,8),
    dims=("time","lat","lon"),
    coords={
        "time": np.arange(10),
        "lat": np.linspace(-90,90,5),
        "lon": np.linspace(0,360,8)
    },
    name="u_wind"
)

DataSet

ds = xr.Dataset(
    {
        "temp": (("lat","lon"), np.random.rand(3,4)),
        "precip": (("lat","lon"), np.random.rand(3,4))
    },
    coords={
        "lat": [10,20,30],
        "lon": [100,110,120,130]
    }
)

print(ds)

大規模データ並列処理

daskで複数ファイルにまたがる大きなファイルを取り扱える。

ds = xr.open_mfdataset("*.nc", chunks={"time":12})