チリ地震津波で観測された津波の高さを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