LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@jp_yen(切手 一円)

BME280 で環境測定 (octave で読み込んでグラフ化)

Last updated at Posted at 2019-03-06

ログを octave でグラフに変換しよう。

ログファイルの 1カラム目が time_t なので、time_t を日付に変換する方法を探す。

T = strftime("%F %T",localtime(1551846602))

C と互換の方法ですね。でもこいつは

octave
octave:25> A = [
>     1551843601, 17.28;
>     1551844201, 17.99;
>     1551844802, 18.32;
>     1551845401, 18.23; ];
octave:26> 
octave:26> T = strftime("%F %T",localtime( A(:,1) ))
T = 2019-03-06 12:40:01
octave:27> T = strftime("%F %T",localtime( A(1,1) ))
T = 2019-03-06 12:40:01
octave:28> T = strftime("%F %T",localtime( A(2,1) ))
T = 2019-03-06 12:50:01

とまぁ、1つめの要素しか読んでくれない。

その他に何か無いか？と探すと datenum() という関数が。
日付と時刻を DateNumber シリアル日付値に変換する関数。

octave
octave:42> datenum('2019-03-03 10:10:10')
ans =     737487.423726852
octave:43> datenum('2019-03-04 10:10:10')
ans =     737488.423726852

どうも一日で 1カウントアップするらしいので time_t に合わせるには 24h * 60m * 60s = 86400 を掛ければ良さそう。
カウントの開始の日付を確認 time_t = 0 の日付を確認して、その日付のシリアル日付値 (DateNumber) を確認。

tcsh
Ras1:~> date --date=@0 +"%F %T"
1970-01-01 09:00:00
octave
octave:49> datenum('1970-01-01 09:00:00')
ans =           719529.375

time_t = 0 の時、DateNumber = 719529.375
関係式としては $ \frac{\mbox{time_t}}{24\times60\times60} = \mbox{datenum()} - 719529.375 $

ということは、日付の形式に直さなくても octave のシリアル日付値に変換すれば OK？

octave
octave:95> A = [1551843601, 17.28;
>     1551844201, 17.99;
>     1551844802, 18.32;
>     1551845401, 18.23 ];
octave:96> 
octave:96> A (:,1) = A(:,1) ./ (24*60*60) + 719529.375
A =

      737490.527789352                 17.28
      737490.534733796                 17.99
      737490.541689815                 18.32
      737490.548622685                 18.23

わかりやすい日付に変換するには

octave
octave:131> datestr(A(:,1) / (24*60*60) + 719529.375, "'yyyy-mm-dd HH:MM:SS'")
ans =

'2019-03-06 12:40:01'
'2019-03-06 12:50:01'
'2019-03-06 13:00:02'
'2019-03-06 13:10:01'

とかできるんですが、文字列なので行列に入れることができません。そして、シリアル値のほうが扱いやすそうなので DateValue で行くことにします。

cvs ファイルを直接読んでグラフにしてみる。
とりあえず気温だけのグラフを書いてみる。

drow_temp.m
logfile = '/var/log/BME280/201903.csv';
logs=dlmread(logfile, ",");

# time_t / 24/60/60 = datenum() -719529.375
log_date = logs(:,1) / (24*60*60) + 719529.375;
log_temp = logs(:,2);

plot(log_date, log_temp, '-x;Temp;');

ylabel('Temperature')
xlabel('Date')

# X軸の表記を日時にする
datetick('x', "%m/%d %H");

drow_temp.jpg

線を引くとやけにかくかくしてたのでマーカーをつけてみた。
うーん、線は適当にしか引いてくれないのね。これはマーカーだけにしたほうが良さそうですね。
X軸は日付で表示させることに成功。さて、これに気圧とか表示させるには Y 軸をもう一本書く必要があるんだけど、はて？

plotyy ってのが使えるらしい。なるほど。
ってことですっ飛ばしてこんな感じになりました。

write_graph.m
# !/usr/bin/octave-cli -qf
# {

 write_graph.m <input.cvs> <output dir>
 
 input.cvs のファイルを読んで output_dir/ の下に jpeg ファイルを出力する
 出力するファイル名は input.cvs のファイル名の拡張子を jpg に変えたもの

# }

format long
clear

# 引数のチェック
# 入力ファイル、出力ディレクトリの存在を確認
# {
if (nargin != 2)
	printf ("%s <input.cvs> <output dir>\n", program_name);
	return
endif
# }

# 引数に指定が無いときのダミー用ファイル
log_file = '/var/log/BME280/2019-03.csv';
graph_dir = tilde_expand("~pi/BME280/");

if (1 <= nargin)
	log_file = tilde_expand(argv(){1})
endif
if (exist(log_file,"file") == 0)
	error("input file not found!: %s", log_file)
	return
endif

if (2 == nargin)
	graph_dir = tilde_expand(argv(){2})
endif
if (exist(graph_dir,"dir") == 0)
#	mkdir(graph_dir)
	error("output directory not found!: %s", graph_dir)
	return
endif
# make path
[path ,name, ext] = fileparts(log_file);
graph_file = fullfile(graph_dir, strcat(name, '.png'));


# CVS ファイルを読む
logs=dlmread(log_file, ",");

# time_t / 24/60/60 = datenum() -719529.375
log_date = logs(:,1) / (24*60*60) + 719529.375;
log_temp = logs(:,2);
log_pres = logs(:,3);
log_humm = logs(:,4);

g_handle = figure('visible', 'on');

### #{
############################
### スプライン曲線で補間 ###
############################
log_date_div = min(log_date):(1/24/60):max(log_date);
log_temp_pp = spline (log_date, log_temp);
log_pres_pp = spline (log_date, log_pres);
log_humm_pp = spline (log_date, log_humm);

[Yaxis,  g_temp,  g_pres ] = plotyy(
	log_date_div, ppval(log_temp_pp,log_date_div),
	log_date_div, ppval(log_pres_pp,log_date_div))
hold on
[~    ,  g_humm,  ~      ] = plotyy( Yaxis(1),
	log_date_div, ppval(log_humm_pp,log_date_div), 0, 0)
set(g_humm, 'color', [0  0  .3])
set(g_temp, 'color', [.3 0  0 ])
set(g_pres, 'color', [0  .3 0 ])
graphs=[g_temp,g_pres,g_humm]
set(graphs, 'marker','.', "markersize", 3, 'linestyle', 'none')
### #}

############################
###   データをプロット   ###
############################
[Yaxis, g_temp, g_pres ] = plotyy(log_date, log_temp, log_date, log_pres);
grid on
grid minor on
hold on
[~     , g_humm, ~      ] = plotyy(Yaxis(1), log_date, log_humm, 0, 0);
# g_humm = plot(log_date, log_humm, 'g');

title(name,
	'fontsize', 14,
	'fontname', 'Noto Mono',
	'color', [0 0 0])
################
# 線の色, 種類 #
################
graphs=[g_temp,g_pres,g_humm];
set(g_temp, 'color', [0.9 0.3 0.3])
set(g_pres, 'color', [0.3 0.9 0.3])
set(g_humm, 'color', [0.3 0.3 0.9])
set(graphs, 'marker','x', "markersize", 3, 'linestyle', 'none')

legend(graphs, 'Temperature [^\\circC]', 'Pressure [hPa]', 'Relative Humidity [%]', 'location', 'northwest')
legend("boxoff")

################
###   Y 軸   ###
################
ylim(Yaxis(1), [floor(min([log_temp;log_humm])/10)*10, ceil(max([log_temp;log_humm])/10)*10])
ylim(Yaxis(2), [floor(min(log_pres)/10)*10, ceil(max(log_pres)/10)*10])

set(Yaxis(1), 'ylabel', 'Temperature[^\\circC] and R.H.[%]', 'ycolor', [0.9 0.3 0.3])
set(Yaxis(2), 'ylabel', 'Atmospheric pressure[hPa]',         'ycolor', [0.3 0.9 0.3])
set(Yaxis, 'fontsize', 14)

################
###   X 軸   ###
################
# ラベルを表示するのは一つだけにする
set(Yaxis, 'XTickLabel', "");

# 書式を日時にする
datetick();		# 日付の自動書式
datetick('x', "%d"); # 指定書式の場合

xlabel('Date')
# X 軸の範囲を揃える
xlim(Yaxis(1), [floor(min(log_date)), ceil(max(log_date))])
xlim(Yaxis(2), [floor(min(log_date)), ceil(max(log_date))])

# 1時間前 = datenum(0,0,0,1)
# 1日前   = datenum(0,0,1)
# xlim(Yaxis(1), [now()-datenum(0,0,0,1), now()])
# xlim(Yaxis(2), [now()-datenum(0,0,0,1), now()])

saveas(g_handle, graph_file)

引数の対応は後ほどきちんと行うとして、とりあえずグラフの作成できた。
赤：温度[℃]、青：湿度[%]、緑：大気圧[hPa]
明るいｘ印が実際の測定値で、濃い点がスプライン曲線で補間した値。6日の夜中から 7日の朝方の温度、湿度が怪しいですね。補間はしちゃ駄目ですね。
2019-03.png
補間部分をコメントアウトし、ｘ印ではなく、点でプロットしたグラフ。12日の湿度が急激に下がったのは前日の雨の後に窓を開けたからです。
2019-03.png

めも

plot の属性一覧
https://octave.org/doc/interpreter/Line-Properties.html#Line-Properties

ax = plotyy(~,~)
set (ax, 'xscale', 'log','xminortick', 'on');
X軸 を対数メモリ、minortick を on にする。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?