概要
この記事では、Octaveを使って不連続に記録されたデータを日付順にプロットする方法について述べます。(たぶんMATLABでも同様に実装できます。試していませんが…)
きっかけ
大学の課題でOctaveを使っていてつまずいたので、対処方法とともにまとめようと思ったから。
環境
・macOS BigSur 11.1
・GNU Octave 6.1.0
・Visual Studio Code 1.52.1
はじめに
気象庁のホームページから過去の気象観測データをCSV形式でダウンロードすることができます。
・過去の気象データ・ダウンロード―気象庁
ここでは、2019年7月1日〜2021年2月1日の山形市の日平均(気温・湿度・気圧)の観測データをダウンロードし、凡例項目などを消してtxtファイルとして保存しています。
対比のために、以下2つのファイルを作成しました。
・data_cont.txt …期間内全てが連続しているデータ(582日分)
・data_disc.txt …連続したデータから~~(適当に)~~ランダムに削除し、不連続にしたデータ(413日分)
これらのデータの中から気温についてグラフにし画像として保存します。
連続したデータの場合
元データとしてdata_cont.txtを用いて、下記のようなコードを作成しました。
clear all
data = load("data_cont.txt"); #データ読み込み
temp = data(:,4); #気温
days = size(temp)
tt = [1:(days)];
# プロット作成
plot(tt, temp)
# 凡例作成
h=legend({"Temperature"},"Location","bestoutside") #凡例位置はプロットにかぶらない最適な位置で、グラフの外側を指定
set(h,"FontName","Times New Roman","FontSize",11)
set(gca,"FontName","Times New Roman","FontSize",11)
xlabel("Number of days elapsed since July 1, 2019(days)")
ylabel("Temperature(degree Celsius)")
# プロットをpngファイルに出力
print("Yamagata2019-21_cont","-dpng")
このコードを実行すると、下の画像が出力・保存されます。
不連続なデータの場合
元データとしてdata_disc.txtを用います。
連続したデータの場合と同様のコードを作成して実行すると下のような画像が出力・保存されます。
何かがおかしいですね…
「はじめに」で述べたようにdata_disc.txtは582日という期間のうち、413日分のデータが含まれているファイルになります。
したがって、経過日数を413日とする上のグラフは図として正しくないです。
対処方法
Octaveにはdatenum()という関数が用意されています。この関数を用いて日付をシリアル日付に変換することができます。
(参照：日付と時刻をシリアル日付値に変換ーMathWorks 日本)
そこで、以下のように記述すると日付の値をx軸に用いる下準備ができます。
Y = data(:,1); #年
M = data(:,2); #月
D = data(:,3); #日
Days = [datenum(Y,M,D)];
次にベクトルの長さを取得し、x軸用に空ベクトルを定義します。
l = length(Days)
tt = [];
ここで、シリアル日付値で最初のデータの日付を表すとどうなるのか確認しておきます。コマンドラインでOctaveを起動し、datenum(年,月,日)のように入力すると知ることができます。
2019年7月1日について求めてみると
octave:1> datenum(2019,7,1)
ans = 737607
となり、737607であることがわかります。
最後にfor文を用いて、先ほど定義した空ベクトルに経過日数を格納していきます。
for i = 1:l
sd = Days(i);
day = sd - 737607; #シリアル日付値で2019年7月1日は737607
tt(i) = day;
end
最初のデータが記録された日を0とすることで経過日数を求めています。
完成形
このようなコードになりました。
clear all
data = load("data_disc.txt"); #データ読み込み
Y = data(:,1); #年
M = data(:,2); #月
D = data(:,3); #日
temp = data(:,4); #気温
Days = [datenum(Y,M,D)]; #各データの年月日をシリアル日付値に変換してベクトルへ格納
l = length(Days) #ベクトルの長さを変数に代入
tt = []; #空ベクトル定義
# x軸プロット用forループ
for i = 1:l
sd = Days(i);
day = sd - 737607; #シリアル日付値で2019年7月1日は737607
tt(i) = day;
end
# プロット作成
plot(tt, temp)
# 凡例作成
h=legend({"Temperature"},"Location","bestoutside") #凡例位置はプロットにかぶらない最適な位置で、グラフの外側を指定
set(h,"FontName","Times New Roman","FontSize",11)
set(gca,"FontName","Times New Roman","FontSize",11)
xlabel("Number of days elapsed since July 1, 2019(days)")
ylabel("Temperature(degree Celsius)")
# プロットをpngファイルに出力
print("Yamagata2019-21_disc","-dpng")
このコードを実行すると、下の画像が出力・保存されます。
正しい経過日数でグラフが描けました。
終わりに
OctaveのTipsや対処方法などは日本語の記事があまりありません。誰かの役に立てば幸いです。
この記事で用いたコードは全てGithubにアップロードしてあります。参考にしたい方はどうぞ。