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|
（略）

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 では使えませんでした．その他にもやってあげる必要がありますので別途まとめました．

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

• 10分角 → 約20km（18km）
• 5分角 → 約10km（9km）
• 1分角 → 約2km（1.8km）
• 30秒角 → 約1km
• 10秒角 → 約300m
• 1秒角 → 約30m

なので，世界の地形データでよく使われている SRTM30 の30は，30秒角を意味していて，大体1kmの空間解像度があると分かります．

QGISでcsvファイルのポイント情報を描画

QGISでcsvのポイント情報を読み込む手順のメモです．ここでは，地震本部の観測施設一覧の内，検潮・津波観測施設（←このリンクは地震本部のWebGISが立ち上がります）を例にしています．

背景は地理院タイルです．
ファイルの読み込みは，レイヤ→レイヤの追加→デリミティッドテキストレイヤの追加…


エンコーディング，ファイル形式（CSV），無視するヘッダー行数（今は無いので0行に修正），XフィールドとYフィールドのカラムを指定します．


[ok]を押すと，CRS（座標系）の設定画面が出てきますので，データのページに「緯度・経度は世界測地系（JGD2000）」と書いてありますので，地理座標系のJGD2000を選択肢ます．と言いたいところですが，見つけるのが大変なので，フィルターのところに入力して検索します．


すると，こんな感じで表示されます．

QGISインストール

別に今更ではあるが，QGISをOS X Sierraにインストールしたのでメモ．

ダウンロードサイトからOS別のファイルをダウンロードする．
「別途 separate installation of dependency framework (GDAL Complete) が必要です。」という注意書きがありますが，この時点では何も考えなくて構いません．KyngChaos QGIS download pageに行くと，CurrentとLong Term Support を選択できますが，分からなければ，Current で良いでしょう．オプションも，NumPy は既にインストーラーに一緒に入っているようなので，ここでは何も考えずに本体のダウンロードだけで進んでしまいましょう．

dmg を開くと，
1 GDAL Complete
2 NumPy
3 matplotlib
4 Install QGIS
と４つのインストーラーが入っていることが分かります．ここでは，上から順に何も考えずに全部インストールしてしまいます．

インストールは簡単に終わりましたが，2.18.14 がリリースされてるぞと警告されてる．．．今持ってきたのになんでだ？と思ったら，10月27日に出たばっかりの模様．Mac版は未だか？

QGISで地理院タイル

QGISで地理院タイルを背景として表示する手順です．

まずは，TileLayer Pluginを使えるようにする必要があります．
まずは，プラグイン→プラグインの管理とインストール　メニューを出します．

検索窓で，tilelayer と入力すると，TileLayer Pluginが表示されると思いますので，インストールします．

そうすると，Webメニューにタイルレイヤプラグインが追加されます．QGIS 2.18.13をインストールして初回に起動した直後にこれをやったら，失敗しました．原因は分かりませんが，プラグインの読み込み自体が出来ていませんでした．プラグインを一旦アンインストールして，QGISを再起動して，問題なく動きました．



タイル情報は別途設定する必要がありますが，TileLayer Plugin作者のminouraさんのGitHubからダウンロードできます．
https://gist.github.com/minorua/7654132

基本的には，
https://gist.github.com/minorua/7654132#file-gsimaps-tsv
を持ってくれば十分です（GSIMaps.tsv）が，[Download ZIP] ボタンで，２つの tsv ファイルを取得できます．

この時，「明確な保存場所」を自分で決めておいてください．ディレクトリを消したり移動したりしてしまうと，QGISが落ちます．

左下の[設定]ボタンから，ZIPを解凍したディレクトリを指定します．この操作は次回以降省略できます．

表示したいレイヤーを選択して[追加]ボタンを押します．

淡色地図で日本付近を拡大すると，以下のような感じで表示できます．

各点の属性情報は，レイヤパネルにあるレイヤで右クリックすると出て来るメニューの「属性テーブルを開く」を選択すると見ることが出来ます．

ただ，これで表が表示されても，どの点の情報なのか分からなくて使い物になりません．．．各点の属性を見るには，上部にある操作メニューアイコンのうち，地物情報表示ボタン（iに矢印）を押して，対象の点をクリックすると右側に地物情報パネルが出てきます．その中の，（アクション）の下のアイコンをクリックすると，その点の地物情報が表示されます．各点の情報を見るにはこっちの方が良いと思います．