2015/09/17

チリ地震津波で観測された津波の高さをGnuplotで書いてみた

2015年9月17日7時54分頃にチリで発生した地震で観測された津波の高さをgnuplotで書いてみました.

データ (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報でプロットしなおしました(余り変わったように見えませんが点が増えてます).
(追記2)
第7報で再プロットです.
set view ,240
でz軸周りに回転して太平洋からの視点に変更しました.

おまけで,軸ラベルも付けています.
set xlabel 'longitude [^o]'
set ylabel 'latitude [^o]'
set zlabel 'Height [m]' rotate by 90