Octaveで間隔がバラバラのデータをプロットする方法

<!-- Place this tag in your head or just before your close body tag. -->
<script async defer src="https://buttons.github.io/buttons.js"></script>

##概要
この記事では、Octaveを使って不連続に記録されたデータを日付順にプロットする方法について述べます。<font color="gray">(たぶんMATLABでも同様に実装できます。試していませんが…)</font>
##きっかけ
大学の課題でOctaveを使っていてつまずいたので、対処方法とともにまとめようと思ったから。
##環境
・macOS BigSur 11.1
・GNU Octave 6.1.0
・Visual Studio Code 1.52.1
##はじめに
気象庁のホームページから過去の気象観測データをCSV形式でダウンロードすることができます。
・[過去の気象データ・ダウンロード―気象庁](https://www.data.jma.go.jp/gmd/risk/obsdl/index.php)

ここでは、2019年7月1日〜2021年2月1日の山形市の日平均(気温・湿度・気圧)の観測データをダウンロードし、凡例項目などを消してtxtファイルとして保存しています。
対比のために、以下2つのファイルを作成しました。

・[data_cont.txt](https://github.com/PyRadiolarus/for-qiita/blob/master/octave/data_cont.txt) …期間内全てが連続しているデータ(582日分)
・[data_disc.txt](https://github.com/PyRadiolarus/for-qiita/blob/master/octave/data_disc.txt) …連続したデータから~~(適当に)~~ランダムに削除し、不連続にしたデータ(413日分)

これらのデータの中から気温についてグラフにし画像として保存します。
##連続したデータの場合
元データとして```data_cont.txt```を用いて、下記のようなコードを作成しました。

```obj-c:data_cont.m
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")
```
このコードを実行すると、下の画像が出力・保存されます。
![2019/07/01~2021/02/01の山形市の日平均気温推移](https://pbs.twimg.com/media/EtdfmztVoAUmzt9?format=png)
##不連続なデータの場合
元データとして```data_disc.txt```を用います。
連続したデータの場合と同様のコードを作成して実行すると下のような画像が出力・保存されます。
![2019/07/01~2021/02/01の山形市の日平均気温推移_何かが違うVer.](https://pbs.twimg.com/media/EtdftoVUcAAgDo0?format=png)
何かがおかしいですね…
「はじめに」で述べたように```data_disc.txt```は582日という期間のうち、413日分のデータが含まれているファイルになります。
したがって、経過日数を413日とする上のグラフは図として正しくないです。
###対処方法
Octaveには```datenum()```という関数が用意されています。この関数を用いて日付をシリアル日付に変換することができます。
(参照：[日付と時刻をシリアル日付値に変換ーMathWorks 日本](https://jp.mathworks.com/help/matlab/ref/datenum.html))
そこで、以下のように記述すると日付の値をx軸に用いる下準備ができます。

```obj-c:
Y = data(:,1); #年
M = data(:,2); #月
D = data(:,3); #日

Days = [datenum(Y,M,D)];
```
次にベクトルの長さを取得し、x軸用に空ベクトルを定義します。

```obj-c:
l = length(Days)
tt = [];
```
ここで、シリアル日付値で最初のデータの日付を表すとどうなるのか確認しておきます。コマンドラインでOctaveを起動し、```datenum(年,月,日)```のように入力すると知ることができます。

2019年7月1日について求めてみると

```
octave:1> datenum(2019,7,1)
ans = 737607
```
となり、```737607```であることがわかります。

最後にfor文を用いて、先ほど定義した空ベクトルに経過日数を格納していきます。

```obj-c:
for i = 1:l
    sd = Days(i);
    day = sd - 737607; #シリアル日付値で2019年7月1日は737607
    tt(i) = day;
end
```
最初のデータが記録された日を0とすることで経過日数を求めています。
##完成形
このようなコードになりました。

```obj-c:data_disc.m
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")
```
このコードを実行すると、下の画像が出力・保存されます。
![2019/07/01~2021/02/01の山形市の日平均気温推移_不連続Ver.](https://pbs.twimg.com/media/EtdfzMqUUAA_Rgi?format=png)
正しい経過日数でグラフが描けました。

##終わりに
OctaveのTipsや対処方法などは日本語の記事があまりありません。誰かの役に立てば幸いです。

この記事で用いたコードは全てGithubにアップロードしてあります。参考にしたい方はどうぞ。

https://github.com/PyRadiolarus/for-qiita/tree/master/octave