gnuplotでDARTの観測記録をプロットしてみる

インドネシアのスマトラ島（Sumatra）沖で，M7.8の地震が2016-03-02 12:49:48 (UTC)に発生しました．幸い横ずれだったようで津波はほとんど発生していません．

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

gnuplotフォントでエラー

自分が管理しているワケじゃない環境で，gnuplot が font 周りでエラーを出したので調査．

エラーメッセージ：
Could not find/open font when opening font "arial", using internal non-scalable font

scalableなフォントがないので，non-scalable なフォントで代用しますよということなので，見た目がきちゃないフォントになってしまいました．警告なのでグラフは書けます．環境は，CentOS 6.6 / gnuplot 4.2 patchlevel 6です．

対処策：
export GDFONTPATH=/usr/share/fonts/liberation
export GNUPLOT_DEFAULT_GDFONT=LiberationSans-Regular
GNUPLOT_DEFAULT_GDFONTはこのディレクトリにあるやつからお好みで選べば良いと思います．サンセリフ体の標準であるこれで良いと思いますが・・・．

参考（というか答えそのものですが．．．）：
http://www5a.biglobe.ne.jp/~nkgwtty/nn_gnuplot.html

gnuplotで使える色

gnuplotで線の色を変えたい場合，

  plot "data" lt 3

などとして（ltをlcにすれば，色だけを指定）指定します．どんな色が使えるかは，test コマンドでplotしてみると分かります（set terminal 毎に異なります）．

set term png
set output "test.png"
test
unset output

で作成したのが以下のプロットです．

ですが，例えば白の線を書きたいと思っても番号がありません（白じゃ見えないじゃないかと言われるかもしれませんが，込み入った線の上から白で線を書きたくてこの発想に至りました）．そのような時は，色の名前で指定することも出来ます．

  plot "data" w l lc rgb "white"

どんな色の名称を使用できるかというのは，

  show colornames

で確認出来ます．

        There are 111 predefined color names:
  white              #ffffff = 255 255 255
  black              #000000 =   0   0   0
  dark-grey          #a0a0a0 = 160 160 160
  red                #ff0000 = 255   0   0
  web-green          #00c000 =   0 192   0
  web-blue           #0080ff =   0 128 255
  dark-magenta       #c000ff = 192   0 255
  dark-cyan          #00eeee =   0 238 238
  dark-orange        #c04000 = 192  64   0
....

冒頭に書いてある通り，111種類定義されているようです．

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

これを

% convert -delay 100 obs*.png obs.gif

でgifアニメーション化

（東経と西経を別々に作成してmontageで連結しています）

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

gnuplotでカラーマップを書く

gnuplotで３次元データを色で描画する方法．

set size sq
set view map
set pm3d explicit map interpolate 10,10
set surface

set title 'x^2+y^2'
set xlabel 'x'
set ylabel 'y'

set xrange [-10:10]
set yrange [-10:10]

set cbtics 25

set palette define (0 "light-green", 1 "yellow", 2 "red")

set output "palette.png"
set term png enh size 640,640 font ",24"
splot x**2-y**2

unset output

gnuplotのpm3dで色を変更する

pm3dで描画した図の色がどぎついと言われてしまったので変更を試みる．

set pm3d map
set isosamples 50
set size sq
unset key
set palette rgbformulae 36,0,-36
splot [-3:3] [-3:3] x*x+y*y-9

こんな感じで青から，黒，赤と正負の値を描けるようになります．

が，クセモノはこの rgbformulae の数字．

デフォルトのまま描画する（set palette しない）とこんな色．

gnuplot> show palette
        palette is COLOR
        rgb color mapping by rgbformulae are 7,5,15
        figure is POSITIVE
        all color formulae ARE NOT written into output postscript file
        allocating ALL remaining color positions for discrete palette terminals
        Color-Model: RGB
        gamma is 1.5

なので，色を反転させようかなと思って，
set palette rgbformulae 15,5,7
としても，全く思ったようにいかないというか，全然違う色になる．

と思ったら，同じ悩みなページ（どせいけいさんき。）を発見．なるほど，rgbformulaeの値をset palette defineの後に得ようというものみたいですが，このページでも指摘されている通り，set palette define だけで私には十分だった．

set pm3d map
set isosamples 50
set size sq
unset key
set palette define (0 "light-green", 1 "yellow", 2 "red")
splot [-3:3] [-3:3] x*x+y*y-9

gnuplotでpm3dを使って2次元のデータも描画する

Gnuplot を使って等高線を書く時は set contour で出来た．最近のGnuplotでは密度グラフが書けるようになったということで，その時のメモ

set view 0,0
3次元グラフを真上から見る

set pm3d explicit map
explicit は，2次元のラインデータなども同時に描画するため
mapは，Y軸メモリを右ではなく左に持ってくるため

unset ztics
しないと，view 0,0 の時にZ軸のラベルが縮退してしまってなんだか分からなくなる．

set cbrange [min:max]
カラーバーのスケールを決める．xrange などと同様に指定がなければ適当にやってくれる．複数のグラフを書くときは設定する方が吉．

splot 'data' notitle w pm3d, 'world_10m_3d.txt' lt -1 notitle w l

world_10m_3d.txt は，Gnuplot's Tips で紹介されていた世界地図描画用のデータを3次元バージョンに変換したもの．と言っても，splotで読めるようになれば良いだけなので，x y に対して， z を某かの固定値で追加してあげれば良い．

＊ 書いてて判明したけど，2次元データをそのまま using で3次元に見立てて上げれば良いだけみたい．

'world_10m.txt' using 1:2:0 lt -1 notitle w l

でも描画された（何故こうなるのかは不明）

参考： - not so Frequently Asked Questions -