アホみたいに時間を浪費したので、メモを日本語でまとめておきます。
今回は、gx1v6グリッドを1度x1度グリッドに変換しました。
ゴール
海洋モデルPOP独自の座標 (グリーンランドに極がある) を、普通の緯度経度格子座標 (1度x1度) に変換する。
問題点
- CMIP5-PMIP3で提供されている出力が、モデル座標のものしかなかった。
- ググったところ、NCLで変換する方法は、わりとたくさん紹介されていた (例）。だが、筆者はNCLは今まで使ったことがなく、このためだけにゼロからNCLを勉強するのは避けたい。
解説している記事を発見
こちらのページにある方法でうまくいった。
以下、日本語でまとめなおしたが、適宜引用元を参照されたし。
手順
グリッド情報を入手
POP-toolという、POPの特殊なグリッドを描画するためのPythonパッケージがある。
今回は、グリッド情報をnetCDF形式で得るために利用する。
まずは、インストール。
conda install -c conda-forge pop-tools
軸情報をダウンロード。
import xarray as xr
import pop_tools
ds = pop_tools.get_grid('POP_gx1v6', scrip=True)
ds.info()
特に指定しなければ、~/.pop_tools以下にインストールされるよう。
ファイル形式(.ieeer8)をnetCDFに変換。
ds.to_netcdf('~/.pop_tools/inputdata/ocn/pop/gx1v6/grid/gx1v6_grid.nc')
変換先の軸情ファイルを作成
NCOのncremapコマンドを使う。
まずは、ncoをインストール。
conda install -c conda-forge nco
これで、ncremapが使えるようになる。
今回は1度x1度の直交座標系に変換したいので、下記のように指定し、grid_1x1_uni_global.ncというグリッド情報ファイルを作成。
ncremap -g grid_1x1_uni_global.nc -G ttl='Equi-Angular grid 180x360'#latlon=180,360#lat_typ=uni
グリッド変換用の重み付けファイルを作成
引き続き、ncremapコマンドを使用して、gx1v6から1度x1度の直交座標系に変換するための重み付けファイルgx1v6_to_1x1_uni_global.ncを作成。
ncremap -a bilinear -s ~/.pop_tools/inputdata/ocn/pop/gx1v6/grid/gx1v6_grid.nc -g grid_1x1_uni_global.nc -m gx1v6_to_1x1_uni_global.nc
グリッド変換
さきほどNCO (nremap) をインストールした際に、同時にインストールされたコマンド群のなかに、ncksがある。
ncksコマンドを使用して、大元のNetCDFファイル (tos_Oclim_CCSM4_piControl_r1i1p1_180001-230012-clim.nc) から必要な変数 (tos) のみを切り出し、tos_gx1v6.ncとして出力。
ncks -v tos tos_Oclim_CCSM4_piControl_r1i1p1_180001-230012-clim.nc tos_gx1v6.nc
cdo setmissval,-999.0 tos_gx1v6.nc out.nc ; mv out.nc tos_gx1v6.nc # ついでに欠損値の変更 (必須ではない)
ncremapコマンドでグリッド変換！！
ncremap -m gx1v6_to_1x1_uni_global.nc tos_gx1v6.nc tos_180x360.nc
ncdump -hとかncviewとかしてみると、グリッド変換できているのが確認できる。やったね！！
軸が(i,j)となっており、(lon,lat)に直したい。今回は、下記の方法で別のファイルで同じ(lon,lat)座標のものから軸情報をうつすことに成功。ncremapはNCOコマンド群に同梱。
ncrename -d i,lon -d j,lat tos_180x360.nc tmp.nc # 軸の名前を変更。
cdo merge tmp.nc ../sic.nc tmp2.nc # 正しい軸情報を持っているファイルといったんマージ。
ncks -v tos tmp2.nc tmp3.nc # 改めて、必要な変数のみ書き出し。
mv tmp3.nc tos_180x360.nc
おまけ
海陸マップを変更した古気候実験の結果に同じ方法を適用したら、グリッド変換後に海岸線付近に変な値が入りました・・・。本来はgrid情報ファイルの海岸線を変更すべきなのだろうが、ひとまず海水温を地表気温を利用して埋めました。