pythonのxarrayは多次元配列を座標付き扱うライブラリで、netcdfにも対応しており、大気海洋モデルの解析に非常に便利である。しかし、緯度経度座標に特化したライブラリではないので、例えば等緯度経度座標のグリッドでも高緯度ではグリッドが小さくなるということを計算に入れてはくれない。それを補うライブラリの一つがxgcmである。これによって微分・積分操作が非常に楽になる(例)。
Version 0.5より、鉛直座標変換が可能になったと知った。
Vertical coordinate transformation is now a part of #xgcm. 🎉
— Julius Busecke (@JuliusBusecke) September 25, 2020
Its also super fast thanks to @numba_jit and @dask_dev.
Try it now!https://t.co/OPYBjPQBtD#openscience #CDSLab @pangeo_data@xarray_dev #opensourcehttps://t.co/2B1EiPLdlf pic.twitter.com/yTiQkySiHi
そこで海洋モデルのσ座標をz座標に変換してみた。単純な線形補間の例である。
uという緯度・経度・σレベルのxarrayの3次元Data arrayを座標変換する。uは
print(u)
結果:
<xarray.DataArray 'um' (lat: 56, lon: 73, lev: 46)>
dask.array<xarray-<this-array>, shape=(56, 73, 46), dtype=float64, chunksize=(56, 73, 46), chunktype=numpy.ndarray>
Coordinates:
* lat (lat) float32 33.98611 34.01389 34.041668 ... 35.48611 35.51389
* lon (lon) float32 138.0 138.02777 138.05556 ... 139.97223 140.0
* lev (lev) float32 1.0 2.0 3.0 4.0 5.0 6.0 ... 42.0 43.0 44.0 45.0 46.0
lonが経度、latが緯度、levがσ座標である。levは
print(u.lev)
結果:
<xarray.DataArray 'lev' (lev: 46)>
array([ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14.,
15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28.,
29., 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., 40., 41., 42.,
43., 44., 45., 46.], dtype=float32)
Coordinates:
* lev (lev) float32 1.0 2.0 3.0 4.0 5.0 6.0 ... 42.0 43.0 44.0 45.0 46.0
σ座標の番号が入っている。
それぞれのσ座標の深さは3次元配列でzzuというxarrayのData arrayに入っている。
print(zzu)
結果:
<xarray.DataArray 'zzu' (lev: 46, lat: 56, lon: 73)>
[188048 values with dtype=float64]
Coordinates:
* lev (lev) float64 1.0 2.0 3.0 4.0 5.0 6.0 ... 42.0 43.0 44.0 45.0 46.0
* lat (lat) float64 33.99 34.01 34.04 34.07 ... 35.43 35.46 35.49 35.51
* lon (lon) float64 138.0 138.0 138.1 138.1 ... 139.9 139.9 140.0 140.0
Attributes:
name: zz
z座標は負で入っている。例えば10番目(0から数えて)の経度・15番目(0から数えて)の緯度の値は
print(zzu.isel(lon=10,lat=15))
結果
<xarray.DataArray 'zzu' (lev: 46)>
array([ -0.996002, -3.486006, -7.470013, -13.446023, -22.908039,
-35.856061, -50.796086, -65.816626, -78.343768, -87.381504,
-92.431898, -95.406424, -98.380951, -101.355474, -104.330004,
-107.304531, -110.279062, -113.253588, -116.228126, -119.202649,
-122.177172, -125.151695, -128.126218, -131.100752, -134.075283,
-137.049809, -140.024332, -142.998855, -145.973385, -148.947916,
-151.922454, -154.896969, -157.8715 , -160.846015, -163.820553,
-166.795091, -169.769606, -172.744137, -175.71866 , -178.69319 ,
-181.667713, -184.642244, -187.616774, -190.591297, -193.565827,
-196.540327])
Coordinates:
* lev (lev) float64 1.0 2.0 3.0 4.0 5.0 6.0 ... 42.0 43.0 44.0 45.0 46.0
lat float64 34.4
lon float64 138.3
Attributes:
name: zz
uとzzuは同じ座標を共有しているが、lon,lat,levの並びは違っている。そこのところはライブラリがよろしくやってくれるようで、後の操作には影響しない。
ここでxgcmが登場で、X座標がlat、Y座標がlon、Z座標がlevであることを定める。cyclicな座標ではないことも宣言する(XやY座標をcyclicにすることもできる)。座標の値はuから取る。
from xgcm import Grid
grid = Grid(u, coords={'Z': {'center': 'lev'},
'X': {'center': 'lon'},
'Y': {'center': 'lat'}},
periodic=False)
print(grid)
結果:
<xgcm.Grid>
Z Axis (not periodic, boundary=None):
* center lev
X Axis (not periodic, boundary=None):
* center lon
Y Axis (not periodic, boundary=None):
* center lat
保存的な座標変換ではグリッドの境界("outer")も定義する必要があるが(例)、今は単純な線形補間なのでグリッドの位置("center")だけで良い。
座標変換は次のようにすれば良い。uをZの方向に[-50,-100,-20]のz座標に変換する。z座標はzzuに入っており、補間方法はlinearである。あえてひねって[-50,-100,-20]と単調増加でない変換をしてみたが、問題ないようである。
uz = grid.transform(u, 'Z', [-10,-500.,-20], target_data=zzu, method='linear')
uzは
print(uz)
結果:
<xarray.DataArray 'um' (lat: 56, lon: 73, zzu: 3)>
dask.array<interp_1d_linear, shape=(56, 73, 3), dtype=float64, chunksize=(56, 73, 3), chunktype=numpy.ndarray>
Coordinates:
* lat (lat) float32 33.98611 34.01389 34.041668 ... 35.48611 35.51389
* lon (lon) float32 138.0 138.02777 138.05556 ... 139.97223 140.0
* zzu (zzu) float64 -10.0 -500.0 -20.0
lon,lat,zzuの配列に変換できている。