EKBO: 11月 2017

2017/11/29

JAGURSで遠地津波計算

津波の伝播計算を自分がやりたいようにやるには，結局一番大変なのは海底地形データを準備することのようです．一例として，遠地地震（チリ地震等）による津波を計算するため，グローバルな地形データを用意する手順をまとめておきます．

まずは，地形データを入手する必要があります．陸域であれば，SRTMやG-DEMが使えますが，海では，GEBCO (General Bathymetric Chart of the Ocean) を使うことが出来ます．ユーザ登録をする必要はありますが，フリーで使えます（詳細なライセンスはサイトを確認してください．以下でダウンロードした zip にも入っています）．

ダウンロードすべきは，"GEBCO_2014 Grid (30 arc-second interval) の User-defined area or global grid | 2D netCDF" です．これにチェックをして，"Add data to basket" を押して，"View basket" に進んで待ちます．分かりやすいインタフェースとは言い難いですが，リクエストしたデータが用意できるまで少し待つ必要があります．データが準備出来たら，ダウンロードすると zip ファイルを取得できます．解凍すると GEBCO_2014_2D.nc が入っていて，これをJAGURSで使える形に変換します．

注意点は以下の通り

極が入っていると計算できないので，必要な領域（太平洋を中心として日本とチリが含まれる範囲等）の切り出します．
この時，日付変更線（東経１８０度と西経１８０度）をまたぐので，東経で１９０度等のように東経で0度から360度までの連続的な座標に変換することも同時に必要です．
JAGURSは水深が正なので，標高データであるGEBCOのデータを正負反転（-1をかける）させます．
NetCDFの標準（？）である nf ではなく，JAGURSは cf を前提にで作られているので変換します．

cfとnfは，grdreformat のオプションで確認することが出来ました．cfは，GMT3の古いフォーマットのfloatで，nfがGMT4のデフォルトのフォーマットのようです．以下は grdreformatの出力の抜粋．

| ID | DESCRIPTION - GMT 3 netCDF legacy formats |

| cf | GMT netCDF format (float) (deprecated) |

----------------------------------------------------------------

| ID | DESCRIPTION - GMT 4 netCDF standard |

| nf | GMT netCDF format (float) (COARDS-compliant) [DEFAULT] |

いくつかやり方はあるはずですが，gmt4を使うなら，

% grdcut GEBCO_2014_2D.nc -Gtmp.grd -R120/300/-60/60

% grdsample tmp.grd -Gtmp2.grd -I2m

% grdmath tmp2.grd NEG = gebco2min.grd=cf

という感じになります．grdcut で領域の切り出しをしています．東端を300とすることで東経の連続データとして変換します．grdsample でも-Rを指定できますが，120/300 のように指定をしても，意図したようには機能してくれませんでした．（必要に応じて，30秒グリッドのオリジナルデータを grdsample で解像度を変更しています．ここでは， -I2m で，２分角に変更しています．）最後に，正負を逆転（grdmath の NEG）しつつ， =cf でフォーマットをcf に変換しています．
※ 領域を Baba et al. (2017) に合わせました．

出来上がった gebco.grd は，JAGRUS のgridfile.dat で指定する地形データとして利用することが出来ます．ちゃんと確認していませんが，どうやらディレクトリは正しく処理できないようなので（たぶん，delimiterの違い？），実行ディレクトリにシンボリックリンクを作成するのが妥当かと思います．

gmt4は，port install gmt4 でインストールしています．

% makecpt -I -Cglobe -Z > gebco.cpt

% pscoast -U -Rgebco2min.grd -JQ150/18 -V -Dc -W1 -B30g30/15g15 -K > gebco2min.ps

% grdimage gebco.grd -J -Rgebco.grd -O -P -Cgebco.cpt >> gebco2min.ps

したのが以下の感じです（90度回転や切り出しは別途しています）．

参考として，２分角メッシュでdt=5秒で24時間分の計算を，MacBook Airでやってみました．OMP_NUM_THREADS=2 で２並列にしています．参照→JAGURS on Mac．

[Timer Output]
+-------------------------------+----------+------------+
|Timer region |Called |Elapsed |
| | |Time[s] |
+-------------------------------+----------+------------+
|All | 1| 5663.732|
|getpar | 1| 0.001|
|read_grid_info | 1| 0.000|
|read_bathymetry_gmt_grdhdr | 1| 0.017|
|read_bathymetry_gmt_grd | 1| 1.331|
|read_friction_gmt_grd | 1| 0.167|
|initl_wfld | 1| 0.715|
|wet_or_dry | 1| 0.103|
|boundary_rwg | 1| 0.001|
|maxgrd_init_rwg | 1| 0.148|
|maxgrd_v_init_rwg | 1| 0.126|
|tgs_open_ct | 1| 0.000|
|hrise_rwg | 12| 0.394|
|drise_rwg | 12| 1.451|
|dump_gmt_nl_vel | 97| 0.002|
|dump_gmt_nl_hgt | 97| 55.068|
|tstep_grid_vel | 17281| 1265.025|
|- fxy_rwg | 17281| 1264.911|
|- mapgrids | 34562| 0.004|
|- outsea_rwg | 34562| 7.162|
|- recheck_wod | 34562| 0.003|
|maxgrd_v_check_nl | 17281| 2804.151|
|tstep_grid_hgt | 17281| 657.405|
|- hxy_rwg | 17281| 650.099|
|maxgrd_check_nl | 17281| 370.801|
|error_check | 17281| 271.747|
|maxgrd_write_gmt | 1| 0.558|
|maxgrd_v_write_gmt | 1| 0.528|
+-------------------------------+----------+------------+

MacProで4スレッド計算だと以下の感じです．
[Timer Output]
+-------------------------------+----------+------------+
|Timer region |Called |Elapsed |
| | |Time[s] |
+-------------------------------+----------+------------+
|All | 1| 3255.667|
（略）

2017/11/16

gdalでnetCDF

MacPorts の gdal は，netCDF 対応は variants になってるので作り直しました．これを機に，もろもろ variants をつけて対応フォーマットを増やしておきます．

port install gdal +expat +netcdf +geos +hdf4 +hdf5 +xerces

これで，海底地形データ（GEBCO）を読むことが出来るようになります．GEBCO は，Geral Bathmetry Chart of the Oceans の略で，最大30秒角（約1km）の海底地形のグリッドデータのことです．遠地の津波を計算する場合にはとりあえず用いられているようです．

持ってきたのは，Gridded bathmetry data の GEBCO_2014 Grid (30 arc-second interval) の2D netCDFです．不思議なことに，INT16 GeoTIFF (data) を持ってこようとチェックすると，
Selected image size 43200 pixels width by 21600 pixels height (Max size allowed 10800 pixels width by 10800 pixels height)
と怒られてしまいます．フォーマットでピクセル数は変わらないと思うのですが．．． 2D netCDF を持ってくれば良さそうです．

% gdalinfo GEBCO_2014_2D.nc
Driver: netCDF/Network Common Data Format
Files: GEBCO_2014_2D.nc
Size is 43200, 21600
.....

INT16で配布されている GEBCO_2014_2D.nc を JAGURS で使えるように cf (Float32な古いnetCDF形式）にします．

% grdreformat GEBCO_2014_2D.nc GEBCO_2014_2D_cf.nc=cf
（これでちゃんと使えるかの確認はこれから）
　　→ cfに変換するだけでは JAGURS では使えませんでした．その他にもやってあげる必要がありますので別途まとめました．

（余談）秒角と距離の対応：
当然ですが，緯度が変わると経度方向の角度と実距離は変化するので固定ではありませんが，緯度方向は一定になっています．ざっくりですが，赤道（か，少しずれた範囲）を基準に考えて，大凡の実距離に換算できます．