Pico Tech • Gnuplot Memo3D

Radiation Detection / English Documents ... Francais　... 日本語ページインデックス / Survey Links ... 掲示板 / 放射線量 / 線量データ / Gnuplot Memo / Peak Fit With Gnuplot

Ｇｎｕｐｌｏｔで、時系列スペクトルファイルなどの３Ｄ表示する

２０１１年３月１４日、２２時頃の福島県大野モニタリングポストの３Ｄスペクトル

環境（背景）放射線のモニタリングなどでは、時系列の多数のスペクトルファイルが作成され、
それらを観察して色々検討したり確認するには、動画やアニメーションがしばしば使われますが、
３Ｄ表示にすると（特に大きな変化などは）とても分かりやすくなることがあります。

そこで、Ｇｎｕｐｌｏｔの使い方にも徐々に慣れてきたので、既に沢山保存されている
各種のスペクトルファイルから、ほぼ自動的に３Ｄスペクトル画像を作成する方法を
考えてみました。

Ｇｎｕｐｌｏｔをまだインストールしてなかったり、どういうものかご存じない方は、
まず、こちらのメモをご覧ください。
http://pico.dreamhosters.com/GnuplotMemo.html

初めは、テレミノＭＣＡで常に測定し続けているうちのＢＧモニタリングデータを使って色々とテストしました。

理想的には、データが入っているフォルダーを指定すると、後はお任せで適当な３Ｄスペクトルを
作ってくれると良いのですが、Ｇｎｕｐｌｏｔが扱える３Ｄグラフのデータを調べると、
未加工のテレミノＭＣＡのスペクトルをそのまま与えて３Ｄ表示をするのは無理みたいなので、
ＰｙｔｈｏｎでＧｎｕｐｌｏｔ用のデータファイルを作成することにしました。

３Ｄ　データの形式

データは、以下のような３列（又はそれ以上）になっているか、
そういう形にどうにか持っていかないとなりません。

「日時」,「チャンネル、又は、エネルギー」,「カウント、又は、ｃｐｓ」

この内、「日時」の部分は、スペクトルデータの場合ファイル名やヘッダー部分に
含まれているものを抽出したりしてそれを各チャンネルのデータに付け加えれば良いでしょう。

Ｓｐｅファイルの様に、チャンネル番号もエネルギーも無い場合には、
その部分も適当に生成して、くっつけないとならないでしょう。

３Ｄ描画用のコマンド部分

現在使っているのは、こんな感じ。
細かな設定は、必要に応じて変更しています。

いくつか注意する点がありますが、
軽く動く（マウスでグラフを動かして、色々な角度から見たり、拡大、縮小したりする）のを
重視する場合には、（特に初期表示の）表示範囲を小さくしたり、ズームして、
データが余り沢山描画されないようにしたり、色をつけずにＭｅｓｈで表示したり、
最初か最後のデータを重複させることにより、「グリッドデータ」として
扱われずに、（Ｍｅｓｈの横線も表示せず）単に各ファイルのスペクトルの線だけ
表示したりするようにすると、良いでしょう。

見栄え重視の場合は、色が付いていた方が、綺麗で、高級そうに見えたりします。

細かな様子を見るには、色が付いていない方が見やすい場合と、
色つきの方が見やすい場合があり、表示の角度や、陰線処理なども
場合によりけりなので、表示方法を変えて、色々な状態で観察すると良さそうです。
これは、グラフのサイズにも言えて、グラフ表示用の窓を大きくして見れば
分かりやすくなるかと言うと、一概にそうともいえず、縦横比やディスプレーの具合にも
左右されたりするので、やっぱり、あれこれ試してみると良いようです。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 
#!/usr/bin/gnuplot --persist

reset

# データを読み込む際の日時のフォーマットを指定。
# 日時のフォーマットは、他の言語などとほぼ同じ。
set timefmt x "%y%m%d%H"
set timefmt x "%Y%m%d_%H%M"

# 横軸が日時であることを指定
set xdata time

# 横軸の日時の表示フォーマットの指定。　お好みに合わせて変更可能
set format x "%Y/%m/%d %H:%M"
#set format x "%H:%M"

# グリッドを表示
set grid

# ｃｓｖなので、データの区切りはカンマを指定
set datafile separator ","

# 線グラフを使うことを指定
set style data lines

#set encoding "utf8"

# 縦軸（Ｚ軸）は対数表示
set log z
#set xrange[-10:3000]
set zrange [1:100000]
set zrange [1:*]

# 横軸（エネルギー）は２０keVから２ＭｅＶまで
set yrange[20:2000]
set yrange[0:500]
#set zrange[1:220]

# Ｚ軸の接地面の調整
set ticslevel 0.01

# 初期表示角度の設定
set view 60,120

# 陰線にする場合の設定
set hidden3d

# グリッドデータで、色付きにする場合の各種設定
# 色の指定は、好きなものを選んで適当に調整。
set pm3d
set palette defined (-3 "blue", 0 "white", 1 "red")
#set palette rgbformulae 33,13,10
set palette model HSV rgbformulae 3.5,2,2
set logscale cb

# ラベルなどの設定
set xlabel "Time"
set xlabel ""
set ylabel "Enregy keV"
set zlabel "Counts" offset graph 0,0,0.6
set xtics  #offset 2

# 目盛りの表示を前に持ってくる
set tics front

# 実際はチャンネルだけど、分かりやすい様に大よそのエネルギーで表示
set ytics ("500" 100,"1000" 200,"1500" 300,"2000" 400)

set zrange [1:*] 
set cbrange [1:*] 

# 初期表示のセット（出来るだけ軽く動く様に、ズームしてデータの少ない状態）
set xrange ["11031200":"11031100"] reverse
#set xrange [*:*] reverse

#splot fn u 2:1:3 t "" lc rgb "#60d4683a"
#set term pngcairo size 800,480 font "MS UI Gothic,10"

# （必要なら）適当なタイトルを付ける
set title 'TITLE'

# "splot" が３Ｄグラフのコマンド。　このコマンドに続けて、データを
# コマンドファイルに含める場合は、ファイル名の代わりに　"-"　を入れる 
splot "-" u 1:2:3 t "" lc rgb "#60d4683a"

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上の内容が入っているファイル：
http://pico.dreamhosters.com/raddata/spehead.plt

上記のコマンドスクリプトの最後の行を抜いて適当な名前で保存し、
それを対話環境で　load "cmd.plt" などと一度呼び出せば、
後は、　splot "data.csv" u 1:2:3 t "" lc rgb "#60d4683a"
とかやれば、良いですし、データの上に、上記のコマンドを全部貼り付けて、
load "data.plt" とかやっても、良いでしょう。

マウスで３Ｄグラフを動かす

対話環境で表示した場合は、マウスでもって、グラフを回転させたり、
角度を変えたり、表示する日時の範囲を拡大/縮小したり、移動させることが出来ます。

マウスの左ボタンをクリックして動かすと、グラフを回転させられます。

マウスホイールを回すと、Ｙ軸（エネルギー）の表示範囲が移動します。
（初期表示より高い／低いエネルギーの方を見るのに便利）

Ｓｈｉｆｔキーを押しながらマウスホイールを回すと、時間軸（Ｘ軸）の表示範囲が移動します。
ＳｈｉｆｔとＣｏｎｔｒｏｌキーを押しながらマウスホイールを回すと時間軸の拡大、縮小が出来ます。
私は、この二つを良く使います。

キーボードで３Ｄグラフを操作する

マウスホイールでの操作が非常に便利ですが、
「対話環境でsplot コマンド使用した場合」であるなら、
キーボード操作でコマンドを実行させて、表示範囲やズームの変更を
させることが出来ます。

ただし、今のところ、この方法は「コマンドとデータが一体化されているファイル」を
load コマンドで呼び出した場合や、エクスプローラーからダブルクリックして
描画した場合などは、使用できていません。

この方法は、　bind コマンドという、一種のキーマクロ指定のコマンド使います。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 
bind 'y' 'gpy=0.2*(GPVAL_Y_MAX-GPVAL_Y_MIN);set yrange[GPVAL_Y_MIN+gpy:GPVAL_Y_MAX-gpy]; replot'
bind 'Y' 'gpy=0.2*(GPVAL_Y_MAX-GPVAL_Y_MIN);set yrange[GPVAL_Y_MIN-gpy:GPVAL_Y_MAX+gpy]; replot'
bind 'alt-y' 'gpy=0.1*(GPVAL_Y_MAX-GPVAL_Y_MIN);set yrange[GPVAL_Y_MIN+gpy:GPVAL_Y_MAX+gpy]; replot'
bind 'alt-Y' 'gpy=0.1*(GPVAL_Y_MAX-GPVAL_Y_MIN);set yrange[GPVAL_Y_MIN-gpy:GPVAL_Y_MAX-gpy]; replot'

bind 'x' 'gpx=0.2*(GPVAL_X_MAX-GPVAL_X_MIN);set xrange[GPVAL_X_MIN+gpx:GPVAL_X_MAX-gpx]; replot'
bind 'X' 'gpx=0.2*(GPVAL_X_MAX-GPVAL_X_MIN);set xrange[GPVAL_X_MIN-gpx:GPVAL_X_MAX+gpx]; replot'
bind 'alt-x' 'gpx=0.1*(GPVAL_X_MAX-GPVAL_X_MIN);set xrange[GPVAL_X_MIN+gpx:GPVAL_X_MAX+gpx]; replot'
bind 'alt-X' 'gpx=0.1*(GPVAL_X_MAX-GPVAL_X_MIN);set xrange[GPVAL_X_MIN-gpx:GPVAL_X_MAX-gpx]; replot'

bind 'z' 'set zrange[GPVAL_Z_MIN*10:GPVAL_Z_MAX*0.1]; replot'
bind 'Z' 'set zrange[GPVAL_Z_MIN*0.1:GPVAL_Z_MAX*10]; replot'
bind 'alt-z' 'set zrange[GPVAL_Z_MIN*10:GPVAL_Z_MAX*10]; replot'
bind 'alt-Z' 'set zrange[GPVAL_Z_MIN*0.1:GPVAL_Z_MAX*0.1]; replot'

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上記のコマンドを実行すると、ｙとシフトｙでＹ軸（エネルギー）の表示範囲の移動を行い、
Ａｌｔ－ｙ、Ａｌｔ－Ｙで、ズームの変更をします。
Ｘ軸（時間軸）、Ｚ軸（カウントかｃｐｓ）についても、同様です。

Ｇｎｕｐｌｏｔの内部変数である　GPVAL_Y_MAX　と　GPVAL_Y_MIN　に、現在のＹ軸の
表示最大値と最小値が入っているので、その差が表示範囲の幅になります。
この表示幅の２０％をｇｐｙという変数に入れ、Ｙ軸の表示範囲を指定しなおし、
グラフを再描画をする、という操作をキーコードに　結びつける（bind）のが最初の行。
２行目は、似たようなことをやってますが、表示範囲を動かす方向が逆。
３行目と４行目は、表示幅を上下それぞれ１０％づつ広げる、というもの。

Ｘ軸についても、やっていることは、全く同じです。
Ｚ軸は、対数表示なので、２０％とか足したり引いたりする代わりに、
単純に１０倍したり、１０で割ったりしています。

これらは、ＳｔａｃｋＯｖｅｒＦｌｏｗに出ていた例を適当にアレンジしました。

本当は、これが、対話環境以外、そして、コマンド一体型データを　load した時でも
使えると嬉しいのですが。

エクスプローラーでダブルクリックした時も、マウス操作が出来るようにする

Ｇｎｕｐｌｏｔは、ちょっと偏屈なところがありますので、Windows版でも
Ｃｙｇｗｉｎ版でも、対話型環境でプロットしないと、マウスが使えません。
（これは、マウスイベントを監視したり登録して、色々動かすプログラムが
対話型の方にあって、グラフ表示用の窓とセットになっていないから？）

そこで、これを回避するには、Ｇｎｕｐｌｏｔスクリプト（規定値では .plt）に
結び付けられているプログラムを、Windows用に色々アレンジしてある
Ｗｇｎｕｐｌｏｔ．ｅｘｅから、Ｕｎｉｘ系やＣｙｇｗｉｎなどとほぼ共通のＧｎｕｐｌｏｔ．ｅｘｅに変更し、
さらに起動時のパラメーターに　－　（半角のハイフン）を加えてやると、
エクスプローラーでダブルクリックなどで、描画する時も対話環境が
Ｄｏｓ窓で開き、目出度くマウスやマウスホイールで、グリグリとグラフを
動かせるようになります。

具体的には、Windowsキーを押して、Regedit と入力して
ＲｅｇＥｄｉｔ．ｅｘｅを起動し、wgnuplot.exe のシェルオープンのところを
変更します。

あるいは、Ｍｉｎｇｗ版のＧｎｕｐｌｏｔを何も変更しないでインストールした場合、
以下の内容のファイルを実行して、レジストリーを変更します。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

Windows Registry Editor Version 5.00

[HKEY_CLASSES_ROOT\gnuplot\shell\open\command]
@="\"C:\\Program Files (x86)\\gnuplot\\bin\\gnuplot.exe\" -p \"%1\" -"

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

上記のファイル：
http://pico.dreamhosters.com/raddata/gnuplot-mod.reg

強いて難点を挙げるなら、グラフを閉じたら、Ｄｏｓ窓のＧｎｕｐｌｏｔも閉じないとならない点。
まあ、Ｑのキーを２回押してエンターを押すだけで両方閉じますので、苦にはなりませんが。

テレミノＭＣＡのスペクトルをＧｎｕｐｌｏｔで使えるデータ・コマンドファイルにする

Ｐｙｔｈｏｎで書いたテキトーなスクリプトは、これです。
http://pico.dreamhosters.com/raddata/prepspe.py

スペクトルの種類や設定などによって、２，３のバリエーションがあって、
基本的にはスペクトルファイルが入っているディレクトリーを指定し、
場合によっては、最初の日付をしていすると、そのディレクトリーの中の
ファイルを全部まとめたデータ・コマンド一体型ファイルを作るか、
一日毎、あるいは１ヶ月ごとの一体型ファイルを生成し、
同時に、それらを使って自動的に画像ファイルを作るコマンドスクリプトも作ります。

テレミノＭＣＡをご利用の方からのリクエストなどがあれば、
もう少し使いやすい形にまとめます。

このスクリプトで作ったデータファイルと、Ｇｎｕｐｌｏｔで出力した各種の３Ｄスペクトルなどは、
ここに入れてありますので、ご覧ください。
http://pico.dreamhosters.com/raddata/ca/

例：

福島県が公開している、原発事故当時のスペクトルデータを使う

データの置き場をちょろちょろ変更するみたいです・・・
現在は、ここにあります。
http://www.atom-moc.pref.fukushima.jp/old/monitoring/monitoring201103/201103_mpdata.html

こっちは、ファイル構造が単純なので、Ｐｙｔｈｏｎではなく、シェルスクリプトでやりました。
このページの一番上のグラフを作ったのが、このスクリプトです。

スクリプトは、１０分値のデータ用のものと、一時間値用のものがあります。

１０分値用：
http://pico.dreamhosters.com/raddata/speplt.sh

一時間値用：
http://pico.dreamhosters.com/raddata/speplt6.sh

ただ、既にデータ・コマンド一体型ファイルも出来上がっていて、公開していますし、
特にこれらを使う必要はないかもしれませんが、シェルスクリプトだけで、
Ａｗｋも何も使わずに、こういう処理をする（あまり良くないかもしれない）例として
出しておきます。

中でcatを使ってしまっている以外は、純粋なシェルスクリプトです。
（catを使わない方法もありますが）

上記のスクリプトで作ったデータファイルや、３Ｄスペクトル画像は、
以下のディレクトリーに入れてあります。（リンクが間違っていたので修正しました。）
http://pico.dreamhosters.com/raddata/ja/2011_spector_all/

出力済みの原発事故直後の３Ｄスペクトル

以下のページに、出力済みの３Ｄスペクトル画像と、そのデータ付きプロットファイルをまとめてあります。
http://pico.dreamhosters.com/FukushimaDisasterSpectrum3D.html

前処理無し、Ｇｎｕｐｌｏｔ以外の外部コマンド無しで直接３Ｄスペクトルを出す

沢山あるスペクトルデータから、全く下準備無しに、かつＡｗｋやＳｅｄ，Ｐｙｔｈｏｎ、といった
外部コマンドに依存せずに、純粋にＧｎｕｐｌｏｔの中だけで完結して３Ｄスペクトルを
描画できると、WindowsやＭａｃやＵｎｉｘ系のどのプラットフォームでも
無変更で動くスクリプトにしやすいです。

で、結論から言うと、現状では、メッシュや色つきなどの「面表示」は、
「グリッドデータ形式」でなければならず、下準備無しでこれを行うのは
不可能であることが分かりました。

ただし、線表示で良いなら、下準備も外部コマンドも使わずに、
いきなりスペクトルデータから３Ｄスペクトルを描画することが可能です。

方法は簡単で、上記のプロットコマンドファイル：
http://pico.dreamhosters.com/raddata/spehead.plt
の最後の行のプロット文の代わりに以下の様なコマンドを付けるだけです。

 
# 色を付けたい場合は、線データをグリッド形式データに変換する文を使うと簡単に出来ます。
# これは、pm3dとか、他の色関連のコマンドと同様、上記のプロットコマンドの前にしていする必要があります。
# チャンネルピッチに合わせて、最初の数字を指定しますが、あんまり細かくすると、遅くなりますし、
# 小さなグラフだと、そんなに細かくしても違いが見えないので意味がありません。
set dgrid3d 999,2,5


# データのあるディレクトリーを指定
dd = '../Rad/201103_spector_all/201103_spc10_kooriyama/'

# この例では、３月のデータだけ表示するので、この様に変数「ｍ」の値を固定していますが、
# ループで回して、複数の月にまたがるデータを出すことも出来ます。
m = 3

# 日付と、時間から、ファイル名を生成する関数。この場合は、福島県が公開している
# １時間値のスペクトルデータのファイル名を作っています。
fn(m,d,h) = sprintf('%sSPC-11%02d%02d%02d00.txt', dd, m, d, h)

# ３Ｄプロット文の本体。
# まず、最初のＦｏｒループで、日付（１１日から１３日）を指定し、入れ子になった次のＦｏｒループで
# 時間（０時から２３時）を指定し、using　の最初の列（Ｘ軸）はsprintfで日付と時間から日時のデータを作り、
# ２列目（Ｙ軸）は、チャンネル番号をデータファイルの行数からcolumn()関数の特別な例である０を与える形で求め、
# ３列目に、各データファイルの中のカウント数が入ります。
splot for [d=11:13] for [hr=0:23] fn(m,d,hr) using (sprintf('11%02d%02d%02d', m, d, hr)):(column(0)):1 title ''

そうすると、こんな感じの色付きの線グラフ群が３Ｄ表示されます。

幾つかの地点の線量の推移と気象データ

山口県萩市の２０１４年の線量と気象の推移 Hagi Yamaguchi2014
２０１４年７月にスパイクのあった福島県広野町東下地区集会所 Azumakudari Fukushima
２０１４年７月に線量が（一時的に？）急降下した飯舘村役場　蕨平自治公民館 Warabitaira Fukushima
上記の浮島、千鳥、そして川崎の線量と、羽田の気象データを組み合わせたグラフ Ukishima Chidori Tonomachi2014

福島第一原発周辺の情報マップ：　Fukushima Daiichi Map Of Rad Info
福島第一原発周辺（２０から３０ｋｍくらい）の情報マップ：　Fukushima Daiichi Second Map Of Rad Info
福島第一原発周辺（３０から４０ｋｍくらい）の情報マップ：　Fukushima Daiichi Third Map Of Rad Info
福島第一原発周辺（４０から５０ｋｍくらい）の情報マップ：　Fukushima Daiichi Fourth Map Of Rad Info
福島第一原発周辺（５０から６０ｋｍくらい）の情報マップ：　Fukushima Daiichi Fifth Map Of Rad Info
鮫川村の情報マップ: Samekawamura Map Of Radinfo

川内原発周辺の情報マップ：　Satsuma Sendai Map Of Rad Info

宇都宮の情報マップ： Utsunomiya Map Of Radinfo
埼玉周辺の情報マップ：　Saitama Map Of Rad Info
松本周辺の情報マップ：　Matsumoto Map Of Rad Info
新潟周辺の情報マップ：　Niigata Map Of Rad Info
仙台周辺の情報マップ：　Sendai Map Of Rad Info

全アメダス、気象台情報マップ：　The Map Of Amedas

Radiation Detection / English Documents ... Francais　... 日本語ページインデックス / Survey Links ... 掲示板 / 放射線量 / 線量データ