2017/12/12

Emacs利用のためのiTerm2設定

iTerm2内のEmacsでMewを快適に使うためのキーバインディングを書いてましたが,自分でも新しいマシンに同じ設定をするのが面倒なのでまとめ.

⌘-d ⇒ M-d
⌘-x ⇒ M-x
⌘-< ⇒ M-<
⌘-> ⇒ M->
⌘-u ⇒ M-u
⌘-w ⇒ M-w
⌘-% ⇒ M-%

(副作用)
⌘-w ・・・ Windowのcloseが出来なくなります.が,ターミナルの場合は,exitでclose出来るので弊害はない.

(禁忌)
⌘-v ・・・ こいつはペーストなのでやらない
⌘-$ ・・・ スクリーンショットの呼び出しと衝突する

具体的な手順は,iTerm2内のEmacsでMewを使うに書いてあります.

2017/12/06

JAGURSのデータをQGISで表示する

JAGURSで扱うデータは,NetCDFのちょっと古いバージョンですが標準フォーマットなので,GMTだけではなくQGISでも表示できるはずです.ということで,やってみました.

まずは,JAGURS用に変換した標高データ(gebco2min.grd)です.レイヤーの追加で,座標参照システムとして,EPSG:4326(WGS84な緯度経度)を選択します.


それだけで,モノクロで表示されます.ただし,水深が正の値なので海が白くなってしまいます.


色を付けるには,プロパティ→スタイルで,レンダータイプを「単バンド疑似カラー」にします.色は,RdYIBuにしてみました.負の値が赤に割り当てられていますが,今は水深が正の値なので反転はさせないでおきます.


可視化結果が以下の通り.色は,,,これが良いと言うわけではないので(お薦めを見つけたら書き直します)自由に選んで下さい.


ここに,計算結果のSD01.zmax.grdをオーバーレイします.レイヤーの追加で同じくEPSG:4326を指定すれば読み込み完了です.このデータは正のデータだけとりあえず見たいので,0から1の間で白から赤にします.



これだと,既に読み込んである地形データが白で上塗りされてしまってもったいない状態です.プロパティ→透過性で,「画面から値を追加」で陸域を透過設定します.また,追加のデータなし値に0を追加して波が届かなかった場所を背景の地形にします.


結局,最大波高の色をYIOrRdにしています.

見やすいとは言い難いですが,以下のようになりました.

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で使える形に変換します.

注意点は以下の通り

  • 極が入っていると計算できないので,必要な領域(太平洋を中心として日本とチリが含まれる範囲等)の切り出します.
  • この時,日付変更線(東経180度と西経180度)をまたぐので,東経で190度等のように東経で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 で,2分角に変更しています.)最後に,正負を逆転(grdmathNEG)しつつ, =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度回転や切り出しは別途しています).

参考として,2分角メッシュでdt=5秒で24時間分の計算を,MacBook Airでやってみました.OMP_NUM_THREADS=2 で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 では使えませんでした.その他にもやってあげる必要がありますので別途まとめました

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

  • 10分角 → 約20km(18km)
  • 5分角 → 約10km(9km)
  • 1分角 → 約2km(1.8km)
  • 30秒角 → 約1km
  • 10秒角 → 約300m
  • 1秒角 → 約30m
なので,世界の地形データでよく使われている SRTM30 の30は,30秒角を意味していて,大体1kmの空間解像度があると分かります.

2017/11/13

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

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

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


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


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


すると,こんな感じで表示されます.