xarrayを使ったACCの計算がうまくいかなかった

xarrayでACCを計算する

gradsが使えなくなるので、ACCもできるだけpythonで計算するようにしているが、できたらnumpyではなくてxarrayをうまく使って計算したい。ところがやってみたら意外と難しかったという話。

時間軸が原因で計算できない

ACC (anomaly correlation coefficient)を計算する場合、予測と観測データの間で相関係数を計算する。従って予測と観測で時間の軸が揃っている必要がある。書いてみるとあたりまえだけど、モデル出力と観測データでは軸の属性が違うことは結構当たり前であって、同じ設定だと思っていても属性が微妙にずれると処理できない。

前処理

例えば時間軸を合わせるならCDOを使って時間軸の設定をまずはやってしまう。

cdo settaxis,2001-01-01,00:00:00,1month infile.nc outfile.nc

これでとにかく2001年から１ヶ月単位のデータとして設定はできる。ただしこれだとxarrayでは読みこんでくれない。理不尽。

xarrayでの読み込み方

CDOで変更した形式はxarrayが対応していないので、時間軸についてはdecodeしないものとして読み込む。予測と観測で軸の名前が違っていたら合わせてしまう。また、１度読み込めたら軸についてはpandasを使って内部で設定しなおしてしまう。

ds1=xr.open_dataset('temp_ens.nc',decode_times=False)
ds2=xr.open_dataset('temp_obs.nc',decode_times=False)

ds1=ds1.rename({'time_counter':'time'})
ds1=ds1.rename({'deptht':'depth'})

ds1["time"] = pd.date_range(start="2001-01-01",periods=276,freq="MS")
ds2["time"] = pd.date_range(start="2001-01-01",periods=276,freq="MS")

ACCの計算

気候値とアノマリー

ACCを計算するには気候値とアノマリーを求める必要があるので、以下のようにxarrayで計算する。時間軸は合わせてしまっているのでとても楽。以下では3Dの気候値を計算するのは計算負荷が高すぎるので、水深150mの水温を補間で前もって計算している。

temp1=ds1['temp'].interp(depth=150)
temp2=ds2['temp'].interp(depth=150)

clm1=temp1.groupby('time.month').mean('time')
clm2=temp2.groupby('time.month').mean('time')

ano1=temp1.groupby('time.month')-clm1
ano2=temp2.groupby('time.month')-clm2

ACCの計算

numpy的にそのまま計算してしまう。ただし月単位なので、xarrayでgroupbyしてから計算している。

num = (ano1 * ano2).groupby('time.month').mean('time')
den = np.sqrt(
    (ano1**2).groupby('time.month').mean('time') *
    (ano2**2).groupby('time.month').mean('time')
)
acc = num / den