 |
| ポスター |
2020/03/15
2020/01/15
M 6.4 - 8km S of Indios, Puerto Rico on 7 Jan. 2020
津波波形データベース構築の一例として，１月７日にプエルトリコで発生した地震による津波のデータ収集と計算をやってみた．
波源断層モデル
Moment Tensor (USGS)
(17.916, 293.187) 13.5 km, Dip: 43 deg., Strike: 268 deg., Rake: -58 deg.
M0 - S 関係 (cite) から，
(L, W) = (20.0073963517812 km, 10.00369817589060 km)
slip amount: 0.500816743451796 m
DEM作成
GEBCO 2019から，当該領域の切り出し (grdcut) を行った．GEBCO_2019.nc は， grdmath を用いて正負を反転してある．
% grdmath GEBCO/GEBCO_2019/GEBCO_2019.nc NEG = GEBCO_2019.nc
% grdcut -R290/300/15/20 $bathymetry/GEBCO_2019.nc -Ggebco.grd
GEBCO 2019 から切り出した地形データでJAGURSを走らせると，北西の領域端にかかっている島のところで発散が生じた．そのため，周囲３分角（15秒角の12メッシュ）は，水深100 mよりも浅い海と陸地を水深 100 mで固定する修正を grdmask と grdmath で行った．下記は， GMT 6.0.0 (Homebrew on OS X) を使用している（ -I15s は， -I15c に相当）．
% cat<<EOF | grdmask -A -Gmask.grd -I15s -Rgebco.grd -N100/100/NaN
290:3 15:3
290:3 19:57
299:57 19:57
299:57 15:3
EOF
（解説）
grdmask で，標準入力から縁取りの４点を入力して， -N 領域外/線上/内側 の値を指定する． -Rgebco.grd で元のDEMと範囲を合わせて，分解能も指定する（-I15s）．さらに，大円ではなく直線で四角にするため -A オプションを付ける．
-A Suppress connecting geographic points using great circle arcs, i.e., connect by straight lines,
unless m or p is appended to first follow meridian then parallel, or vice versa.
波源断層モデル
Moment Tensor (USGS)
(17.916, 293.187) 13.5 km, Dip: 43 deg., Strike: 268 deg., Rake: -58 deg.
M0 - S 関係 (cite) から，
(L, W) = (20.0073963517812 km, 10.00369817589060 km)
slip amount: 0.500816743451796 m
DEM作成
GEBCO 2019から，当該領域の切り出し (grdcut) を行った．GEBCO_2019.nc は， grdmath を用いて正負を反転してある．
% grdmath GEBCO/GEBCO_2019/GEBCO_2019.nc NEG = GEBCO_2019.nc
% grdcut -R290/300/15/20 $bathymetry/GEBCO_2019.nc -Ggebco.grd
GEBCO 2019 から切り出した地形データでJAGURSを走らせると，北西の領域端にかかっている島のところで発散が生じた．そのため，周囲３分角（15秒角の12メッシュ）は，水深100 mよりも浅い海と陸地を水深 100 mで固定する修正を grdmask と grdmath で行った．下記は， GMT 6.0.0 (Homebrew on OS X) を使用している（ -I15s は， -I15c に相当）．
% cat<<EOF | grdmask -A -Gmask.grd -I15s -Rgebco.grd -N100/100/NaN
290:3 15:3
290:3 19:57
299:57 19:57
299:57 15:3
EOF
% grdmath gebco.grd mask.grd LT mask.grd MUL 0 NAN gebco.grd AND = masked_gebco.grd=cf
（解説）
grdmask で，標準入力から縁取りの４点を入力して， -N 領域外/線上/内側 の値を指定する． -Rgebco.grd で元のDEMと範囲を合わせて，分解能も指定する（-I15s）．さらに，大円ではなく直線で四角にするため -A オプションを付ける．
-A Suppress connecting geographic points using great circle arcs, i.e., connect by straight lines,
unless m or p is appended to first follow meridian then parallel, or vice versa.
験潮所
ITIC (International Tsunami Information Center) のアラートで 0.02 mの津波が観測されたようなので，UNESCO IOC Sea Level Station Monitoring Facility の event からMagueyes Island PR staion のデータを "show data" から取得． 200〜2000秒のバンドパスフィルタをかけた．
(lat, lon) = (17.9701, -67.0464) だったが，DEMに重ね合わせると，ずれていたので， (17.9636, 292.9536) に移動した．
JAGURS計算
設定した断層パラメータで， dt=1.0 s で計算した．
2018/11/09
直交座標系での計算
JAGURSで直交座標系（CARTESIAN）で計算するには，Makefileで
#CARTESIAN=ON
を有効にするために#のコメントアウトを外します．
従って，球座標（JAGURSのデフォルト）で作成したバイナリとは別になりますので時々で使い分ける場合はファイル名を変えるなり，ディレクトリを別にするなりしておく必要があります．
OpenTSUNAMI （杞憂プロジェクト）はぢめます．詳しくは後日．
#CARTESIAN=ON
を有効にするために#のコメントアウトを外します．
従って，球座標（JAGURSのデフォルト）で作成したバイナリとは別になりますので時々で使い分ける場合はファイル名を変えるなり，ディレクトリを別にするなりしておく必要があります．
OpenTSUNAMI （杞憂プロジェクト）はぢめます．詳しくは後日．
Open TSUNAMI（オープン津波）とか称してオープンソースソフトであるJAGURSとオープンデータ（広義）を使って色々やってみることにしました． https://t.co/hsjagQXTkg— 山本直孝 (@naotakayamamoto) 2017年12月21日
2018/10/27
ガウス分布初期波源計算
JAGURSでは，Gaussian分布を海底地殻変動（初期水位）として計算することが出来ます．が，，，マニュアルには球座標使用の場合しか書かれていないので少し補足です．
tsun.par に
init_disp_gaussian =1
とすることで使えます．で，gaussianというファイル名でパラメータを用意します．が，単位に注意が必要です．
球座標の場合：
&gaussian h0=波高 lon_o=経度[度] lat_o=緯度[度] L=半値幅[km] /
直交座標の場合：
&gaussian h0=波高 lon_o=X座標 lat_o=Y座標 L=半値幅[座標系の単位] /
のように，m（メートル）が単位になっているUTM等を利用している場合，Lで指定する半値幅はmが単位になります．一方で，球座標の場合，座標は「度」単位，半値幅はkm単位で固定です．
なお，波高に負を指定することで凹んだ初期水位を表現することも可能です．実装は mod_init_disp_gaussian.f90 に入っていますので必要な関数に書き直しちゃうのも簡単かなと．
JAGURS V0407以前では，直交座標の場合には波高（h0）が指定できない（1 mに固定）ので注意が必要です．
tsun.par に
init_disp_gaussian =1
とすることで使えます．で，gaussianというファイル名でパラメータを用意します．が，単位に注意が必要です．
球座標の場合：
&gaussian h0=波高 lon_o=経度[度] lat_o=緯度[度] L=半値幅[km] /
直交座標の場合：
&gaussian h0=波高 lon_o=X座標 lat_o=Y座標 L=半値幅[座標系の単位] /
なお，波高に負を指定することで凹んだ初期水位を表現することも可能です．実装は mod_init_disp_gaussian.f90 に入っていますので必要な関数に書き直しちゃうのも簡単かなと．
JAGURS V0407以前では，直交座標の場合には波高（h0）が指定できない（1 mに固定）ので注意が必要です．
2018/01/09
地形モデル
遠地津波計算をするためにGEBCOのデータを変換する手順は以前書きましたので，今回は近地津波の計算をするための直交座標系での地形データの作成手順を書いておきます．
元データは内閣府の南海トラフの巨大地震モデル検討会で作成された地形データがG空間情報センターからダウンロード出来ます．無償にはなっていますが，ダウンロードするためにはアカウントを作成する必要があります．地形データ以外の粗度や堤防，初期水位データ等もダウンロード出来ます．．
例として，JAGURS で使える形式に変換してみます．
手順は大きく分けで３つ．(1)まずはデータをダウンロードして，(2)まともなGISフォーマットに変換して，(3)最後に必要なフォーマットへの変換をします．
地形データは，ここからダウンロードできます．利用規約はもちろん確認が必要ですが，「計算入力データ（地形，粗度，堤防データ）の説明」を確認して下さい．ここに何系がどの領域に対応しているか書かれています．例えば四国を計算する場合は第04系になるのでその地形データをダウンロードします．ダウンロードしたデータにはメタ情報が何も書かれていないので，別途「計算範囲設定」を確認する必要があります．下記の ncols, nrows, xllcorner, yllcorner の値は，エクセルのメッシュ個数（X方向，Y方向），南西端の位置（X座標，Y座標）に対応します．
データは 10f8.2 で記述してあるということなので，Fortranで読み込むのは簡単ですが， f8.2 の場合，標高（陸域）が1000mを超える場合に，符号も含めると -1000.00 となり8桁になってしまうため，前の数字とスペースなしで連結されてしまいます．Fortranなら読めるでしょって言ってしまえばそれまでなのですが，一般的なアスキーフォーマットと考えた場合，スペースで区切られていないと使えません．というわけで， 10f8.2 で read して， 10f9.2 で write します．さらに，アスキーフォーマットとしてGISデータになるようにヘッダを付与します．プログラム（reformat.f）も掲載しておきます．
ncols 720
nrows 540
xllcorner -730200
yllcorner -857700
cellsize 2430
NODATA_value -9999
ヘッダの中身は，さほど難しくなく，ncolsはX方向のメッシュ数，nrowsはY方向のメッシュ数，xllcorner は南西端のX座標，yllcorner は南西端のY座標，cellsize は解像度（上記例は 2430m），NODATA_value は不要ですがお決まりの -9999 を設定しています．
reformat.f:
PROGRAM reformat
real*8 data(10)
character(len=128) :: fname
read(*,'(a)') fname
open(10,file=fname)
do
read(10,'(10f8.2)',end=100) data
write(*,'(10f9.2)') data
end do
100 END
作成したアスキーファイル（bathy.SD01.asc）をJAGURS用のcf形式に変換します．
日本海講側は，中央防災会議のものが公開されています（古いです）．
元データは内閣府の南海トラフの巨大地震モデル検討会で作成された地形データがG空間情報センターからダウンロード出来ます．無償にはなっていますが，ダウンロードするためにはアカウントを作成する必要があります．地形データ以外の粗度や堤防，初期水位データ等もダウンロード出来ます．．
例として，JAGURS で使える形式に変換してみます．
手順は大きく分けで３つ．(1)まずはデータをダウンロードして，(2)まともなGISフォーマットに変換して，(3)最後に必要なフォーマットへの変換をします．
地形データは，ここからダウンロードできます．利用規約はもちろん確認が必要ですが，「計算入力データ（地形，粗度，堤防データ）の説明」を確認して下さい．ここに何系がどの領域に対応しているか書かれています．例えば四国を計算する場合は第04系になるのでその地形データをダウンロードします．ダウンロードしたデータにはメタ情報が何も書かれていないので，別途「計算範囲設定」を確認する必要があります．下記の ncols, nrows, xllcorner, yllcorner の値は，エクセルのメッシュ個数（X方向，Y方向），南西端の位置（X座標，Y座標）に対応します．
データは 10f8.2 で記述してあるということなので，Fortranで読み込むのは簡単ですが， f8.2 の場合，標高（陸域）が1000mを超える場合に，符号も含めると -1000.00 となり8桁になってしまうため，前の数字とスペースなしで連結されてしまいます．Fortranなら読めるでしょって言ってしまえばそれまでなのですが，一般的なアスキーフォーマットと考えた場合，スペースで区切られていないと使えません．というわけで， 10f8.2 で read して， 10f9.2 で write します．さらに，アスキーフォーマットとしてGISデータになるようにヘッダを付与します．プログラム（reformat.f）も掲載しておきます．
ncols 720
nrows 540
xllcorner -730200
yllcorner -857700
cellsize 2430
NODATA_value -9999
ヘッダの中身は，さほど難しくなく，ncolsはX方向のメッシュ数，nrowsはY方向のメッシュ数，xllcorner は南西端のX座標，yllcorner は南西端のY座標，cellsize は解像度（上記例は 2430m），NODATA_value は不要ですがお決まりの -9999 を設定しています．
reformat.f:
PROGRAM reformat
real*8 data(10)
character(len=128) :: fname
read(*,'(a)') fname
open(10,file=fname)
do
read(10,'(10f8.2)',end=100) data
write(*,'(10f9.2)') data
end do
100 END
作成したアスキーファイル（bathy.SD01.asc）をJAGURS用のcf形式に変換します．
% grdreformat bathy.SD01.asc=gd bathy.SD01.grd=cf
これで，ネスティングしない場合は計算できるようになりました．ネスティングする場合は，メッシュの細かい調整が必要なので後日掲載します（一応，手元では動くデータは作れましたがこれで良いか確認します）．日本海講側は，中央防災会議のものが公開されています（古いです）．
2017/12/20
JAGURS on Mac (High Sierra)
*** 最近のMacPorts はproj がバージョン5になって，ディレクトリ構造が変わったので Makefile も更新しています．最新の記事を参照して下さい． ***
JAGURSをMacで使う手順は，JAGURS on Mac でまとめましたが，手元がHigh Sierraになったので仕切り直し．
ライブラリ関係は，MacPorts で入れてます．
netcdf-fortran の variants を +gcc7 にしました．
HPC for Mac OS X は，High Sierra 版がないので，Sierra 対応な 7.1 を入れました．大事なことを書き忘れていた気がします．JAGURSは，V0407 です．Makefile は前回の記事から変更していません．
JAGURSをMacで使う手順は，JAGURS on Mac でまとめましたが，手元がHigh Sierraになったので仕切り直し．
ライブラリ関係は，MacPorts で入れてます．
- NetCDF
- netcdf @4.4.1.1_1+dap+netcdf4 (active)
- sudo port install netcdf +dap +netcdf4
- netcdf-fortran @4.4.4_1+gcc7 (active)
- sudo port install netcdf-fortran +gcc7
- PROJ4
- proj @4.9.3_0 (active)
- sudo port install proj
- FFTW3
- fftw-3 @3.3.5_1+gfortran (active)
netcdf-fortran の variants を +gcc7 にしました．
HPC for Mac OS X は，High Sierra 版がないので，Sierra 対応な 7.1 を入れました．大事なことを書き忘れていた気がします．JAGURSは，V0407 です．Makefile は前回の記事から変更していません．
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にしています．
見やすいとは言い難いですが，以下のようになりました．
まずは，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で使える形に変換します．
注意点は以下の通り
という感じになります．grdcut で領域の切り出しをしています．東端を300とすることで東経の連続データとして変換します．grdsample でも-Rを指定できますが，120/300 のように指定をしても，意図したようには機能してくれませんでした．（必要に応じて，30秒グリッドのオリジナルデータを grdsample で解像度を変更しています．ここでは， -I2m で，２分角に変更しています．）最後に，正負を逆転（grdmath の NEG）しつつ， =cf でフォーマットをcf に変換しています．
※ 領域を Baba et al. (2017) に合わせました．
まずは，地形データを入手する必要があります．陸域であれば，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|
（略）
[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 では使えませんでした．その他にもやってあげる必要がありますので別途まとめました．
（余談）秒角と距離の対応：
当然ですが，緯度が変わると経度方向の角度と実距離は変化するので固定ではありませんが，緯度方向は一定になっています．ざっくりですが，赤道（か，少しずれた範囲）を基準に考えて，大凡の実距離に換算できます．
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を選択肢ます．と言いたいところですが，見つけるのが大変なので，フィルターのところに入力して検索します．
すると，こんな感じで表示されます．
背景は地理院タイルです．
ファイルの読み込みは，レイヤ→レイヤの追加→デリミティッドテキストレイヤの追加…
エンコーディング，ファイル形式（CSV），無視するヘッダー行数（今は無いので0行に修正），XフィールドとYフィールドのカラムを指定します．
[ok]を押すと，CRS（座標系）の設定画面が出てきますので，データのページに「緯度・経度は世界測地系（JGD2000）」と書いてありますので，地理座標系のJGD2000を選択肢ます．と言いたいところですが，見つけるのが大変なので，フィルターのところに入力して検索します．
すると，こんな感じで表示されます．
2017/09/09
GMT on Mac
JAGURS のファイルが GMT グリッドファイルということなので，いよいよちゃんと GMT を使うべく手元にインストール．MacPorts にあるようなのでそのまま使います．新しいの入れておけば良いだろうと思って， gmt5 をインストールしたけど， grdinfo コマンドが無い．．．
% port contents gmt5 | grep grdinfo
/opt/local/share/doc/gmt/html/_sources/grdinfo.rst.txt
/opt/local/share/doc/gmt/html/grdinfo.html
/opt/local/share/doc/gmt/man/man1/grdinfo.1.gz
じゃぁ，と思って gmt4 をインストールしたら，変わったディレクトリにコマンドが入ってました．
% port contents gmt4 | grep grdinfo
/opt/local/lib/gmt4/bin/grdinfo
/opt/local/share/doc/gmt4/html/man/grdinfo.html
/opt/local/share/doc/gmt4/man/man1/grdinfo.1
なので，最後の仕上げに
export PATH=$PATH:/opt/local/lib/gmt4/bin
を追加して完成．
第一歩として試したかった，JAGURS 付属の海底地形データの情報を取得．zmax を使って，計算ステップ間隔を決定します．
% grdinfo bathy.SD01.grd
bathy.SD01.grd: Title:
bathy.SD01.grd: Command: grdsample -V -R140/147:30/34/43 -I18c /home/babat/data/bathymetry_thk/bathy.thk.27s.grd -Gbathy.whole.grd=cf
bathy.SD01.grd: Remark:
bathy.SD01.grd: Gridline node registration used
bathy.SD01.grd: Grid file format: cf (# 10) GMT netCDF format (float) (deprecated)
bathy.SD01.grd: x_min: 140 x_max: 147.5 x_inc: 0.005 name: x nx: 1501
bathy.SD01.grd: y_min: 34 y_max: 43 y_inc: 0.005 name: latitude [degrees_north] ny: 1801
bathy.SD01.grd: z_min: -2096.24316406 z_max: 9788.53125 name: z
bathy.SD01.grd: scale_factor: 1 add_offset: 0
% port contents gmt5 | grep grdinfo
/opt/local/share/doc/gmt/html/_sources/grdinfo.rst.txt
/opt/local/share/doc/gmt/html/grdinfo.html
/opt/local/share/doc/gmt/man/man1/grdinfo.1.gz
じゃぁ，と思って gmt4 をインストールしたら，変わったディレクトリにコマンドが入ってました．
% port contents gmt4 | grep grdinfo
/opt/local/lib/gmt4/bin/grdinfo
/opt/local/share/doc/gmt4/html/man/grdinfo.html
/opt/local/share/doc/gmt4/man/man1/grdinfo.1
なので，最後の仕上げに
export PATH=$PATH:/opt/local/lib/gmt4/bin
を追加して完成．
第一歩として試したかった，JAGURS 付属の海底地形データの情報を取得．zmax を使って，計算ステップ間隔を決定します．
% grdinfo bathy.SD01.grd
bathy.SD01.grd: Title:
bathy.SD01.grd: Command: grdsample -V -R140/147:30/34/43 -I18c /home/babat/data/bathymetry_thk/bathy.thk.27s.grd -Gbathy.whole.grd=cf
bathy.SD01.grd: Remark:
bathy.SD01.grd: Gridline node registration used
bathy.SD01.grd: Grid file format: cf (# 10) GMT netCDF format (float) (deprecated)
bathy.SD01.grd: x_min: 140 x_max: 147.5 x_inc: 0.005 name: x nx: 1501
bathy.SD01.grd: y_min: 34 y_max: 43 y_inc: 0.005 name: latitude [degrees_north] ny: 1801
bathy.SD01.grd: z_min: -2096.24316406 z_max: 9788.53125 name: z
bathy.SD01.grd: scale_factor: 1 add_offset: 0
2017/08/12
JAGURS on Mac
JAGURS（高性能津波計算コード）の作者の馬場さんから，Macで動いたなら手順を公開して欲しいという依頼がありましたので，まとめておきます．
（シリアル版かつ手抜きなので順次更新します）
（追記：シリアル版といいつつ，OpenMPは有効に出来るようでした．というわけで，MacOSな１台のマシンでやるにはOpenMPで十分．）
【試した環境】
【依存ライブラリ】
NetCDF，PROJ4，FFTW3が必要ということなので，MacPorts でお手軽に入れました．なので，インストール先は全て /opt/local になりますので，それに合わせてMakefileを修正して下さい．NetCDFはFortranから使えるようにするため，netcdf-fortranパッケージもインストールして下さい．かつ，FFTW3もFortranから使えるようにするため，gfortran variant などを付けておいて下さい．
（というのを考えるのが面倒なので，コンパイルする手順も別途作成します）
% port install szip
しました．
【コンパイル】
まずは，シリアル版．Makefile.SC_ICEをベースにするようにマニュアルに書いてあるけど，Makefileと同じ内容なので，Makefileでスタートする．
【罠】
フォートランのプログラムファイルが .f90 で提供されています．一方，gnu系（？）のコンパイラでは，.f90の場合，preprocessorが動いてくれません！
#ifdef MULTI
1
Warning: Illegal preprocessor directive
とか言われます．これは， .F90 にすることで解決するそうです．まさか．．．
というわけで，ファイル名を .F90 にしてあげて（厳密には，preprocessorを使うソースだけで良いですが面倒なので全部），Makefileの %.o: %.f90 を %.o: %.F90 にします．
（gfortran に -cpp オプションを渡してあげれば良いだけでした）
【Makefile修正箇所】
FC=/usr/local/bin/gfortran
PROJ4_DIR=/opt/local
CC=/usr/local/bin/gcc
CFLAGS=-g -I${PROJ4_DIR}/include
BASE=-cpp
FFTW3_INCLUDE_DIR=/opt/local/include
FFTW3_LIB=-lfftw3
OPT=-O2 -fopenmp
NETCDF=/opt/local
LIBS=-L$(NETCDF)/lib -lnetcdff -lnetcdf -L/opt/local/lib -lhdf5_hl -lhdf5 -lcurl -lsz -lproj $(FFTW3_LIB)
# MPI=ON
修正済Makefileはここに置いておきます．
＊追加注（2017/10/04）
% make
/usr/local/bin/gcc -c -g -I/opt/local/include mapproject.c
FATAL:/opt/local/bin/../libexec/as/x86_64/as: I don't understand 'm' flag!
make: *** [mapproject.o] Error 1
HPC compilerの /usr/local/bin/gcc を使う場合，as が /opt/local/bin/ になってしまうと，この不具合に当たるようです．PATHの順番を /usr/bin の方が先になっている必要があります．
（シリアル版かつ手抜きなので順次更新します）
（追記：シリアル版といいつつ，OpenMPは有効に出来るようでした．というわけで，MacOSな１台のマシンでやるにはOpenMPで十分．）
【試した環境】
- Mac OS 10.12 on Mac pro
- MacPorts + xcode
- HPC Mac OS X compiler
【依存ライブラリ】
NetCDF，PROJ4，FFTW3が必要ということなので，MacPorts でお手軽に入れました．なので，インストール先は全て /opt/local になりますので，それに合わせてMakefileを修正して下さい．NetCDFはFortranから使えるようにするため，netcdf-fortranパッケージもインストールして下さい．かつ，FFTW3もFortranから使えるようにするため，gfortran variant などを付けておいて下さい．
（というのを考えるのが面倒なので，コンパイルする手順も別途作成します）
- NetCDF
- netcdf @4.4.1.1_0+dap+netcdf4 (active)
- sudo port install netcdf +dap +netcdf4
- netcdf-fortran @4.4.4_0+gcc6 (active)
- sudo port install netcdf-fortran +gcc6
- (High) Sierraは HPC compilerがgcc7なので
+gcc7 かも - PROJ4
- proj @4.9.3_0 (active)
- sudo port install proj
- FFTW3
- fftw-3 @3.3.5_0+gfortran (active)
- sudo port install fftw-3 +gfortran 2018.11.12加筆（書き漏れ）
% port install szip
しました．
【コンパイル】
まずは，シリアル版．Makefile.SC_ICEをベースにするようにマニュアルに書いてあるけど，Makefileと同じ内容なので，Makefileでスタートする．
（gfortran に -cpp オプションを渡してあげれば良いだけでした）
【Makefile修正箇所】
- gfortran で .F90 なファイル（.f90と区別される）にpreprocessorを適用するには， -cpp オプションを指定する．
- gfortran で OpenMP を使うには， -fopenmp オプションを指定する．
FC=/usr/local/bin/gfortran
PROJ4_DIR=/opt/local
CC=/usr/local/bin/gcc
CFLAGS=-g -I${PROJ4_DIR}/include
BASE=-cpp
FFTW3_INCLUDE_DIR=/opt/local/include
FFTW3_LIB=-lfftw3
OPT=-O2 -fopenmp
NETCDF=/opt/local
LIBS=-L$(NETCDF)/lib -lnetcdff -lnetcdf -L/opt/local/lib -lhdf5_hl -lhdf5 -lcurl -lsz -lproj $(FFTW3_LIB)
# MPI=ON
修正済Makefileはここに置いておきます．
＊追加注（2017/10/04）
% make
/usr/local/bin/gcc -c -g -I/opt/local/include mapproject.c
FATAL:/opt/local/bin/../libexec/as/x86_64/as: I don't understand 'm' flag!
make: *** [mapproject.o] Error 1
HPC compilerの /usr/local/bin/gcc を使う場合，as が /opt/local/bin/ になってしまうと，この不具合に当たるようです．PATHの順番を /usr/bin の方が先になっている必要があります．
2016/07/30
Mariana諸島の地震による津波
マグニチュード： 7.7
発生日時： 2118 UTC JUL 29 2016 （7月30日 6時18分 日本時間）
位置： 北緯18.5，東経145.8（Mariana islands; マリアナ諸島）
USGS
発生日時： 2118 UTC JUL 29 2016 （7月30日 6時18分 日本時間）
位置： 北緯18.5，東経145.8（Mariana islands; マリアナ諸島）
USGS
No #tsunami threat from 7.7 magnitude #earthquake in the Mariana Islands (too deep to pose a threat) #PTWC https://t.co/kIUg9dSMaT— NWS PTWC (@NWS_PTWC) 2016年7月29日
津波の危険はありませんが，DARTでは観測されていました．
2016/05/03
2016/4/29メキシコ沖の地震による津波
2016/4/29 01:33 (UT) にメキシコ沖で発生した地震による津波のデータを見てみる．
NAOTIDEを使って潮汐補正をした結果が以下の図です．
震央に近いDART 43413の波形
NAOTIDEを使って潮汐補正をした結果が以下の図です．
NAOTIDEを用いて潮汐補正した波形
津波が観測されている付近を拡大
2016/05/02
NAOTIDEJを使った潮汐計算
津波の観測データ（波高や海底水圧）を見ると，最初に見えてくるのは海洋潮汐による１日に２周期の波です．
これを取り除く方法はいくつかあるようですが，海洋潮汐モデルを用いて必要な地点の予測計算をして引き算する方法があります．ここでは，国立天文台水沢で公開されているNAOTIDEJの短周期海洋潮汐モデル（NAO99Jb；日本周辺モデル）を使った潮汐補正を試してみました．日本周辺モデルは以下の範囲と解像度です．
対応範囲：110°E-165°E, 20°N-65°N
解像度：5 arcmin
NAOTIDEJのREADMEより：
DARTのように海水面上にあるGPSブイのように，海面高度（波高）を計測している場合は，itmode=1を使うようです．S-netやDONET等海底水圧計の場合は，itmode=2 を使うようです．
プログラムと海洋潮汐モデルは以下のファイルをダウンロードしました．
海洋潮汐予測プログラム：naotidej990909.tar.gz
海洋潮汐モデル：nao99Jb.tar.gz + nao99L.tar.gz
海洋潮汐モデル：nao99Jb_gc.tar.gz + nao99L_gc.tar.gz
注） READMEに従ってテスト計算するには，_gc が必要です．
naotestj.f が付属しているのでそれをベースにして，観測点の位置と
x = 146.192d0 ! East longitude in degree
y = 40.302d0 ! North latitude in degree
計算する時刻
c Start epoch
iyear1 = 2016 ! year
imon1 = 4 ! month
iday1 = 17 ! day
ihour1 = 0 ! hour
imin1 = 0 ! minute
c
c End epoch
iyear2 = 2016 ! year
imon2 = 4 ! month
iday2 = 19 ! day
ihour2 = 0 ! hour
imin2 = 0 ! minute
，出力間隔（DARTに合わせて15分）
c Output data interval in minute
dt = 15.d0 ! in minute
を修正して実行しました．NAOTIDEJの出力は，（何処かの）平均海水面からの潮位分なのでDARTの水深から直接引き算（naotestj.f はcm単位の出力で，DARTはm単位なので注意）すれば，潮汐補正されたデータを得ることが出来ます．
注） 時刻はDARTもNAOTIDEJもUTなのでそのまま比較できますが，JSTなデータを扱う場合は注意．
観測（紫）と計算（水色）（縦軸のスケールは合わせてあります．１ticが0.2m）を比べると計算したNAOTIDEJの方が少し小さい振幅なので，結果（緑）に少し潮汐が残っています．
気象庁が設置したのDART 21346のデータのプロット
これを取り除く方法はいくつかあるようですが，海洋潮汐モデルを用いて必要な地点の予測計算をして引き算する方法があります．ここでは，国立天文台水沢で公開されているNAOTIDEJの短周期海洋潮汐モデル（NAO99Jb；日本周辺モデル）を使った潮汐補正を試してみました．日本周辺モデルは以下の範囲と解像度です．
対応範囲：110°E-165°E, 20°N-65°N
解像度：5 arcmin
日本周辺モデルの対応範囲
NAOTIDEJのREADMEより：
- itmode = 1 : geocentric tideを計算します。geocentric tideは海洋潮汐と荷重潮汐の和です。海面高度計データの補正の場合はこれを使用します。
- itmode = 2 : 海底を基準とした純粋な海洋潮汐を計算します。海底圧力計データ等の補正の場合はこれを使用します。
DARTのように海水面上にあるGPSブイのように，海面高度（波高）を計測している場合は，itmode=1を使うようです．S-netやDONET等海底水圧計の場合は，itmode=2 を使うようです．
プログラムと海洋潮汐モデルは以下のファイルをダウンロードしました．
海洋潮汐予測プログラム：naotidej990909.tar.gz
海洋潮汐モデル：nao99Jb.tar.gz + nao99L.tar.gz
海洋潮汐モデル：nao99Jb_gc.tar.gz + nao99L_gc.tar.gz
注） READMEに従ってテスト計算するには，_gc が必要です．
naotestj.f が付属しているのでそれをベースにして，観測点の位置と
x = 146.192d0 ! East longitude in degree
y = 40.302d0 ! North latitude in degree
計算する時刻
c Start epoch
iyear1 = 2016 ! year
imon1 = 4 ! month
iday1 = 17 ! day
ihour1 = 0 ! hour
imin1 = 0 ! minute
c
c End epoch
iyear2 = 2016 ! year
imon2 = 4 ! month
iday2 = 19 ! day
ihour2 = 0 ! hour
imin2 = 0 ! minute
，出力間隔（DARTに合わせて15分）
c Output data interval in minute
dt = 15.d0 ! in minute
を修正して実行しました．NAOTIDEJの出力は，（何処かの）平均海水面からの潮位分なのでDARTの水深から直接引き算（naotestj.f はcm単位の出力で，DARTはm単位なので注意）すれば，潮汐補正されたデータを得ることが出来ます．
注） 時刻はDARTもNAOTIDEJもUTなのでそのまま比較できますが，JSTなデータを扱う場合は注意．
DARTの観測データ（紫；左軸）から，NAOTIDEJで計算した潮汐（水色；右軸）を引き算した結果（緑；右軸）
観測（紫）と計算（水色）（縦軸のスケールは合わせてあります．１ticが0.2m）を比べると計算したNAOTIDEJの方が少し小さい振幅なので，結果（緑）に少し潮汐が残っています．
潮汐を引き算した結果のみのプロット
2016/03/03
gnuplotでDARTの観測記録をプロットしてみる
インドネシアのスマトラ島（Sumatra）沖で，M7.8の地震が2016-03-02 12:49:48 (UTC)に発生しました．幸い横ずれだったようで津波はほとんど発生していません．
一方で，NDBC DART（Deep-ocean Assessment and Reporting of Tsunamis）では４点でトリガがかかっているのでその記録を gnuplot で波形を書いてみました．
DART観測点 23401 のデータを使ってみます．今回のイベントのデータはこのリンクから持って来ることが出来ます．データが記載されているtextboxから単にコピペして， DART23401.dat として保存しました．
フォーマットは，データの先頭に以下のように書いてあります．
#YY MM DD hh mm ss T HEIGHT
#yr mo dy hr mn s - m
2016 03 02 16 59 00 2 3469.152
.....
========== DART23401.gp ===========
set size ratio 0.4 # 時刻（x軸）に長くする．
# データを読み込む時の時刻のフォーマット．timefmtに使える記号は help timefmt で
# 確認できます．
set xdata time
set timefmt "%Y %m %d %H %M %S"
# x軸の目盛の書き方も timefmt と同じやり方で指定します．\n で改行も出来ます．
set xlabel 'Time [UT]'
set format x "%m/%d\n%H:%M"
set ylabel 'meters'
set output "DART23401.png"
set term png enh size 1080,480 # size ratio 0.4にしたので横長の画像サイズ
plot "DART23401.dat" u 1:8 lt 0 w l notitle \
, "< ../DART/readData.pl DART23401.dat" in 0 u 1:8 lc 6 lt 5 t '15-min', \
"" in 1 u 1:8 lc 2 lt 4 t '1-min' \
, "" in 2 u 1:8 lc 7 lt 5 t '15-sec'
unset output
# readData.pl は，
# T = Measurement Type;
# where 1 = 15-minute measurement;
# 2 = 1-minute measurement;
# 3 = 15-second measurement
# の３つに分割して（間に空行を２つ入れる）出力するためのスクリプトです．
# 色分けしないなら， plot "DART23401" u 1:8 lt 0 w l だけでプロット出来ます．
x軸（時刻）を指定してプロットしたい場合は， set xrange を使いますが，フォーマットは timefmt で設定した書式に従って，xrange を設定します．例えば，地震動等の影響で見えている最初の振動を拡大するには以下のようにします．この書式は，format x のものではなく，読み込むときに設定する timefmt です．
set xrange ["2016 03 02 13 04 00":"2016 03 02 13 08 00"]
他の３点のデータも下記からダウンロード出来ます．
DART 56003
DART 23227
DART 56001
全データを一定のxrangeにして縦軸のスケールも合わせて，長周期成分を取り除いたデータを axis x1y2 で追加して作成したpng を montage で合成しました．
montage DART?????.png -geometry 540x240 DART.png
<< help timefmt より抜粋>>
Format Explanation
%d day of the month, 1--31
%m month of the year, 1--12
%y year, 0--99
%Y year, 4-digit
%j day of the year, 1--365
%H hour, 0--24
%M minute, 0--60
%s seconds since the Unix epoch (1970-01-01, 00:00 UTC)
%S second, integer 0--60 on output, (double) on input
%b three-character abbreviation of the name of the month
%B name of the month
No tsunami observed. Warnings issued by Indonesia and Marine Warnings by Australia (for Christmas and Cocos Islands)....
Posted by US NWS Pacific Tsunami Warning Center on 2016年3月2日
一方で，NDBC DART（Deep-ocean Assessment and Reporting of Tsunamis）では４点でトリガがかかっているのでその記録を gnuplot で波形を書いてみました．
DART観測点 23401 のデータを使ってみます．今回のイベントのデータはこのリンクから持って来ることが出来ます．データが記載されているtextboxから単にコピペして， DART23401.dat として保存しました．
フォーマットは，データの先頭に以下のように書いてあります．
#YY MM DD hh mm ss T HEIGHT
#yr mo dy hr mn s - m
2016 03 02 16 59 00 2 3469.152
.....
========== DART23401.gp ===========
set size ratio 0.4 # 時刻（x軸）に長くする．
# データを読み込む時の時刻のフォーマット．timefmtに使える記号は help timefmt で
# 確認できます．
set xdata time
set timefmt "%Y %m %d %H %M %S"
# x軸の目盛の書き方も timefmt と同じやり方で指定します．\n で改行も出来ます．
set xlabel 'Time [UT]'
set format x "%m/%d\n%H:%M"
set ylabel 'meters'
set output "DART23401.png"
set term png enh size 1080,480 # size ratio 0.4にしたので横長の画像サイズ
plot "DART23401.dat" u 1:8 lt 0 w l notitle \
, "< ../DART/readData.pl DART23401.dat" in 0 u 1:8 lc 6 lt 5 t '15-min', \
"" in 1 u 1:8 lc 2 lt 4 t '1-min' \
, "" in 2 u 1:8 lc 7 lt 5 t '15-sec'
unset output
# readData.pl は，
# T = Measurement Type;
# where 1 = 15-minute measurement;
# 2 = 1-minute measurement;
# 3 = 15-second measurement
# の３つに分割して（間に空行を２つ入れる）出力するためのスクリプトです．
# 色分けしないなら， plot "DART23401" u 1:8 lt 0 w l だけでプロット出来ます．
x軸（時刻）を指定してプロットしたい場合は， set xrange を使いますが，フォーマットは timefmt で設定した書式に従って，xrange を設定します．例えば，地震動等の影響で見えている最初の振動を拡大するには以下のようにします．この書式は，format x のものではなく，読み込むときに設定する timefmt です．
set xrange ["2016 03 02 13 04 00":"2016 03 02 13 08 00"]
他の３点のデータも下記からダウンロード出来ます．
DART 56003
DART 23227
DART 56001
全データを一定のxrangeにして縦軸のスケールも合わせて，長周期成分を取り除いたデータを axis x1y2 で追加して作成したpng を montage で合成しました．
montage DART?????.png -geometry 540x240 DART.png
<< help timefmt より抜粋>>
Format Explanation
%d day of the month, 1--31
%m month of the year, 1--12
%y year, 0--99
%Y year, 4-digit
%j day of the year, 1--365
%H hour, 0--24
%M minute, 0--60
%s seconds since the Unix epoch (1970-01-01, 00:00 UTC)
%S second, integer 0--60 on output, (double) on input
%b three-character abbreviation of the name of the month
%B name of the month
2015/10/07
7万3000年前の津波
約7万3000年前の火山噴火による山体崩壊によって巨大な津波が発生していたというニュースがありました．場所は，アフリカ西海岸沖のカーボベルデ諸島のフォゴ島のフォゴ山（Pico do Fogo）．遡上高が800フィート（244メートル）というのは凄いですね．
今でも噴火を繰り返していて，火山衛星画像データベース（産総研）でも観測されています．ASTER高温領域検出システムと火山衛星画像データベースのKMLからFogo山に視点を持ってきたKMLを作成しました．
ちなみに，火山噴火による津波は日本でも起きていて，200年以上前の雲仙岳の噴火では対岸の肥後・天草に50mとも言われる津波が発生していました．
今でも噴火を繰り返していて，火山衛星画像データベース（産総研）でも観測されています．ASTER高温領域検出システムと火山衛星画像データベースのKMLからFogo山に視点を持ってきたKMLを作成しました．
ちなみに，火山噴火による津波は日本でも起きていて，200年以上前の雲仙岳の噴火では対岸の肥後・天草に50mとも言われる津波が発生していました．
2015/09/18
チリ津波で観測された津波高をGnuplotで書いてみた（２）
ちょいちょい修正したので書き直し．その１の続きです．
set zrange [0:5]
set xrange [-180:-20]
set yrange [-65:45]
set ticslevel 0
# 任意の場所に任意の文字（ここでは時刻）を表示
set label ' 1416 UTC THU SEP 17 2015' at -180,45,5
set grid
set xlabel 'longitude [^o]'
set ylabel 'latitude [^o]'
set zlabel 'Height [m]' rotate by 90
set view ,330
# unset border
# mxticsの数字は，xticsの間にいくつ点をうつかという意味（間隔ではない）
set xtics 30 out
set mxtics 3
set ytics 30 out
set mytics 3
set ztics border
h1(x)=(x<1)? x:0 # 1mより低い場合を抽出
h3(x)=(x<3)? x:0 # 3mより低い場合を抽出
set term png enh font "Times,18" size 720,720
set output "obs18.png"
# ./obs.pl obs18.dat はデータ整形スクリプトで lon lat height の数字を出力
splot "world_10m.txt" u 1:2:(0) w lines lt -1 lw 0 notitle \
, "< ./obs.pl obs18.dat" u 1:2:3 w impulse lt 7 lw 3 t '3m<=' \
, "" u 1:2:(h3($3)) w impulse lt 4 lw 3 t '<3m ' \
, "" u 1:2:(h1($3)) w impulse lt 6 lw 3 t '<1m '
unset output
これを
（東経と西経を別々に作成してmontageで連結しています）
set zrange [0:5]
set xrange [-180:-20]
set yrange [-65:45]
set ticslevel 0
# 任意の場所に任意の文字（ここでは時刻）を表示
set label ' 1416 UTC THU SEP 17 2015' at -180,45,5
set grid
set xlabel 'longitude [^o]'
set ylabel 'latitude [^o]'
set zlabel 'Height [m]' rotate by 90
set view ,330
# unset border
# mxticsの数字は，xticsの間にいくつ点をうつかという意味（間隔ではない）
set xtics 30 out
set mxtics 3
set ytics 30 out
set mytics 3
set ztics border
h1(x)=(x<1)? x:0 # 1mより低い場合を抽出
h3(x)=(x<3)? x:0 # 3mより低い場合を抽出
set term png enh font "Times,18" size 720,720
set output "obs18.png"
# ./obs.pl obs18.dat はデータ整形スクリプトで lon lat height の数字を出力
splot "world_10m.txt" u 1:2:(0) w lines lt -1 lw 0 notitle \
, "< ./obs.pl obs18.dat" u 1:2:3 w impulse lt 7 lw 3 t '3m<=' \
, "" u 1:2:(h3($3)) w impulse lt 4 lw 3 t '<3m ' \
, "" u 1:2:(h1($3)) w impulse lt 6 lw 3 t '<1m '
unset output
これを
% convert -delay 100 obs*.png obs.gif
でgifアニメーション化
2015/09/17
チリ地震津波で観測された津波の高さをGnuplotで書いてみた
2015年9月17日7時54分頃にチリで発生した地震で観測された津波の高さをgnuplotで書いてみました．
（追記1）
第6報でプロットしなおしました（余り変わったように見えませんが点が増えてます）．
（追記2）
第7報で再プロットです．
set view ,240
でz軸周りに回転して太平洋からの視点に変更しました．
おまけで，軸ラベルも付けています．
set xlabel 'longitude [^o]'
set ylabel 'latitude [^o]'
set zlabel 'Height [m]' rotate by 90
データ (obs2.dat)
# lat lon m
-26.4 -70.6 1.09
-25.4 -70.5 0.25
-36.7 -73.1 0.42
-20.5 -73.4 0.06
-26.3 -80.1 0.50
-33.0 -71.6 1.62
-33.6 -78.8 0.97
-30.0 -71.3 3.11
-34.6 -72.0 0.46
-27.1 -70.8 0.50
-33.6 -71.6 0.86
-26.7 -74.0 0.12
# CHANARAL CL 26.4S 70.6W 0043 1.09M/ 3.6FT 30
# TALTAL CL 25.4S 70.5W 0018 0.25M/ 0.8FT 38
# TALCAHUANO CL 36.7S 73.1W 0020 0.42M/ 1.4FT 84
# DART 32401 20.5S 73.4W 0020 0.06M/ 0.2FT 48
# SAN FELIX CL 26.3S 80.1W 0019 0.50M/ 1.6FT 06
# VALPARAISO CL 33.0S 71.6W 2350 1.62M/ 5.3FT 40
# JUAN FERNANDEZ 33.6S 78.8W 0008 0.97M/ 3.2FT 08
# COQUIMBO CL 30.0S 71.3W 2339 3.11M/10.2FT 20
# BUCALEMU CL 34.6S 72.0W 2334 0.46M/ 1.5FT 24
# CALDERA CL 27.1S 70.8W 2343 0.50M/ 1.6FT 30
# SAN ANTONIO CL 33.6S 71.6W 2358 0.86M/ 2.8FT 18
# DART 32402 26.7S 74.0W 2326 0.12M/ 0.4FT 36
world_10m.txt は，昔の記事でも使わせて頂いたGnuplot Tips集経由でダウンロードしたものを利用しています．スクリプトは，同じくGnuplot Tips集の3次元棒グラフのページを参考（いや，コピペして修正しただけです）にさせて頂きました．ただし，borderを消しちゃうとどこだか分からなくなっちゃうので，水平面(座標が書かれる面)の位置の変更を参考にticslevelを0にしてxyplaneに地図(z=0)が入るようにしています．
set zrange [0:*]
set xrange [-110:-20]
set yrange [-55:35]
set ticslevel 0
# unset border
set xtics 30
set ytics 30
set ztics border
h1(x)=(x<1)? x:0 # 1mより低い場合を抽出
h3(x)=(x<3)? x:0 # 3mより低い場合を抽出
set term png enh font "Times,18" size 720,720
set output "obs1.png"
splot \
"world_10m.txt" using 1:2:(0) w lines lt -1 lw 0 notitle \
, "obs2.dat" using 2:1:3 w impulse lt 7 lw 3 t '3m<=' \
, "" using 2:1:(h3($3)) w impulse lt 4 lw 3 t '<3m ' \
, "" using 2:1:(h1($3)) w impulse lt 6 lw 3 t '<1m '
unset output
上記のスクリプトで書いたのがこれです．
（補足）
splotで3次元プロットをする場合，xy平面（xyplane）はz軸の適度なところに割り当てられるそうです．それは，以下の式で定められるそうです．
ticslevel = (xyplaneのz座標 - z軸の最小値) / (z軸の最小値 - z軸の最大値)
ここで，xyplaneのz座標が設定したい位置ですね．ticslevelはデフォルトでは0.5になっているので，z軸の最小値より下になります．z軸の最小値にしたい場合は ticslevel=0で出来るようです．
（追記1）
第6報でプロットしなおしました（余り変わったように見えませんが点が増えてます）．
第7報で再プロットです．
set view ,240
でz軸周りに回転して太平洋からの視点に変更しました．
おまけで，軸ラベルも付けています．
set xlabel 'longitude [^o]'
set ylabel 'latitude [^o]'
set zlabel 'Height [m]' rotate by 90
登録:
投稿 (Atom)