xgcmを使ったσ座標からz座標への鉛直座標変換

pythonの[xarray](https://xarray.pydata.org/en/stable/)は多次元配列を座標付き扱うライブラリで、netcdfにも対応しており、大気海洋モデルの解析に非常に便利である。しかし、緯度経度座標に特化したライブラリではないので、例えば等緯度経度座標のグリッドでも高緯度ではグリッドが小さくなるということを計算に入れてはくれない。それを補うライブラリの一つが[xgcm](https://xgcm.readthedocs.io/en/latest/)である。これによって微分・積分操作が非常に楽になる([例](https://xgcm.readthedocs.io/en/latest/example_mitgcm.html))。

[Version 0.5](https://xgcm.readthedocs.io/en/latest/whats-new.html)より、鉛直座標変換が可能になったと知った。
<blockquote class="twitter-tweet"><p lang="en" dir="ltr">Vertical coordinate transformation is now a part of <a href="https://twitter.com/hashtag/xgcm?src=hash&amp;ref_src=twsrc%5Etfw">#xgcm</a>. 🎉<br><br>Its also super fast thanks to <a href="https://twitter.com/numba_jit?ref_src=twsrc%5Etfw">@numba_jit</a> and <a href="https://twitter.com/dask_dev?ref_src=twsrc%5Etfw">@dask_dev</a>.<br><br>Try it now!<a href="https://t.co/OPYBjPQBtD">https://t.co/OPYBjPQBtD</a><a href="https://twitter.com/hashtag/openscience?src=hash&amp;ref_src=twsrc%5Etfw">#openscience</a> <a href="https://twitter.com/hashtag/CDSLab?src=hash&amp;ref_src=twsrc%5Etfw">#CDSLab</a> <a href="https://twitter.com/pangeo_data?ref_src=twsrc%5Etfw">@pangeo_data</a><a href="https://twitter.com/xarray_dev?ref_src=twsrc%5Etfw">@xarray_dev</a> <a href="https://twitter.com/hashtag/opensource?src=hash&amp;ref_src=twsrc%5Etfw">#opensource</a><a href="https://t.co/2B1EiPLdlf">https://t.co/2B1EiPLdlf</a> <a href="https://t.co/yTiQkySiHi">pic.twitter.com/yTiQkySiHi</a></p>&mdash; Julius Busecke (@JuliusBusecke) <a href="https://twitter.com/JuliusBusecke/status/1309540690209255426?ref_src=twsrc%5Etfw">September 25, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>


そこで海洋モデルのσ座標をz座標に変換してみた。単純な線形補間の例である。

uという緯度・経度・σレベルのxarrayの3次元Data arrayを座標変換する。uは

```python
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は

```python
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に入っている。

```python
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から数えて)の緯度の値は

```python
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から取る。

```python
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")も定義する必要があるが([例](https://xgcm.readthedocs.io/en/latest/transform.html))、今は単純な線形補間なのでグリッドの位置("center")だけで良い。

座標変換は次のようにすれば良い。uをZの方向に[-50,-100,-20]のz座標に変換する。z座標はzzuに入っており、補間方法はlinearである。あえてひねって[-50,-100,-20]と単調増加でない変換をしてみたが、問題ないようである。

```code
uz = grid.transform(u, 'Z', [-10,-500.,-20], target_data=zzu, method='linear')
```

uzは

```python
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の配列に変換できている。