気象予報値取得APIを利用して機械学習による電気使用量を予測する手順をまとめてみた

# はじめに

Web上で継続的に公開している電力使用量予測について、新型コロナウイルス感染症の影響からか予測の精度がイマイチになってきた感じがしているところです。

[電力使用量予測 predicted by blueOmega](https://blueomega.jp/ppc_plot.html)

そろそろ再学習させてみようかなと思ったところで、先日開発した気象予報値取得APIを利用する方法を試してみたのでその手順をまとめてみます。
ちなみに、気象予報値取得APIの紹介記事は以下をご参照下さい。

[気象予報値などをAPIで取得できるサービスのテスト版を公開してみた](https://qiita.com/mix_dvd/items/c3c40f8ed6445a16a3fe)

# 学習用データの取得

[でんき予報 | 中国電力ネットワーク](https://www.energia.co.jp/nw/jukyuu/)

上記ページの「過去の電気使用実績」から2019年、2020年の実績をダウンロードします。

[気象庁 | 過去の気象データ・ダウンロード](https://www.data.jma.go.jp/gmd/risk/obsdl/index.php)

上記ページより、広島と松江における１時間ごとの気温データを2019年4月から2020年12月の期間分取得します。

１年ごとにファイルを分けて、以下のファイル名で保存します。

- data-hiroshima-2019.csv
- data-hiroshima-2020.csv
- data-matsue-2019.csv
- data-matsue-2020.csv

# ライブラリ読込

以下のライブラリを使用しました。

```py:
# データ編集用
import math
from glob import glob
import pandas as pd
import requests
import io
from datetime import datetime as dt
from datetime import timedelta

# 機械学習用
from sklearn import model_selection
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.neural_network import MLPRegressor
```

# 関数の宣言

## 気温の読込

```py:
def read_temp(city):
    
    df_tmp = pd.DataFrame()
    
    files = glob("data-{}-*.csv".format(city))
    files.sort()
    
    for f in files:
        df_tmp = pd.concat([df_tmp, pd.read_csv(f,encoding="Shift_JIS",skiprows=4)])
    
    df_tmp = df_tmp.reset_index(drop=True)
    df_tmp.columns = ["DATETIME","TEMP","品質情報","均質番号"]
    
    df_tmp.DATETIME = df_tmp.DATETIME.map(lambda _: pd.to_datetime(_))
    df_tmp.index = df_tmp.DATETIME
    
    return df_tmp
```

# 特徴量の追加

以前公開した以下の記事の機能を実装するために関数を作成します。

https://qiita.com/mix_dvd/items/01230e19bedab8e53ccc

```py:
def make_radian_row(pca_result):
    rad = []
    for r in pca_result:
        rad.append(math.degrees(math.atan(r[0]/r[1])))

    return rad
```

# データ読込

```py:
# 電力データ読込
files = glob("juyo-*.csv")
files.sort()

df_kw = pd.DataFrame()

for f in files:
    df_kw = pd.concat([df_kw, pd.read_csv(f,encoding="Shift_JIS",skiprows=2)])

df_kw = df_kw.reset_index(drop=True)
df_kw["MW"] = df_kw["実績(万kW)"] * 10

df_kw.index = df_kw.index.map(lambda _: pd.to_datetime(df_kw.DATE[_] + " " + df_kw.TIME[_]))

# 気温データ読込

# 広島
df_tmp_hiroshima = read_temp("hiroshima")

# 松江
df_tmp_matsue = read_temp("matsue")
```

# 学習データ作成

```py:
# データの複製
df = df_kw.copy()

# 広島の気温を結合
df["TMP_hiroshima"] = df_tmp_hiroshima.TEMP
# 松江の気温を結合
df["TMP_matsue"] = df_tmp_matsue.TEMP

# 月、週、時間の値を個別に取得
df["MONTH"] = df.index.month
df["DAY"] = df.index.day
df["WEEK"] = df.index.weekday
df["HOUR"] = df.index.hour

# 使用するデータを指定
cols = ['MW', 'TMP_hiroshima', 'TMP_matsue', 'MONTH', 'DAY', 'WEEK', 'HOUR']
df = df[cols]

# null値の削除
df = df.dropna()

# one-hotエンコーディング
cols = ["MONTH","WEEK","HOUR"]
for col in cols:
    df = df.join(pd.get_dummies(df[col], prefix=col))
```

## 学習データと検証データに分割

```py:
# 2020年12月以降のデータを検証用として取得
df_valid = df["2020-12-01":]
df = df[:"2020-11-30"]
```

## 説明変数と目的変数の取得

```py:
# 説明変数に使用するデータを指定
x_cols = ['TMP_hiroshima', 'TMP_matsue', 'DAY']

# 月
for i in range(12):
    x_cols.append("MONTH_{}".format(i+1))

# 曜日
for i in range(7):
    x_cols.append("WEEK_{}".format(i))

# 時間
for i in range(24):
    x_cols.append("HOUR_{}".format(i))

X = df[x_cols]
y = df["MW"]
```

# 学習

```py:
# 学習とテスト用に分割
X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, test_size=.2, random_state=42)

# 特徴量を追加
pca = PCA(n_components=2)
pca.fit(X_train)

X_train["rad"] = make_radian_row(pca.transform(X_train))
X_test["rad"] = make_radian_row(pca.transform(X_test))

# モデルの作成と学習
model = MLPRegressor(hidden_layer_sizes=(45,100,30), random_state=42)
model.fit(X_train, y_train)

print(model.score(X_test, y_test))
```

実行結果

```
0.8989258758760301
```

# 検証

学習データに含まれないテストデータで精度を確認。

```py:
X_valid = df_valid[x_cols]
y_valid = df_valid["MW"]

X_valid["rad"] = make_radian_row(pca.transform(X_valid))

X_valid = scaler.transform(X_valid)

print(model.score(X_valid, y_valid))
```

実行結果

```
0.5613505868156264
```

低っ！

可視化してみる。
![Unknown.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/114240/317cc595-5067-39b1-231b-04d3c081c2bb.png)
最初の方はあってそうだけど、12月の月末から全然だめですね...ま、しかたない。

# 予測

気象予報値を取得するため、以下のサイトのAPIを利用。
http://iot.blueomega.jp/sensingplaza_next/

```py:
x_cols = ['TMP_hiroshima', 'TMP_matsue', 'DAY', 'MONTH_1', 'MONTH_2', 'MONTH_3', 'MONTH_4', 'MONTH_5', 'MONTH_6', 'MONTH_7', 'MONTH_8', 'MONTH_9', 'MONTH_10', 'MONTH_11', 'MONTH_12', 'WEEK_0', 'WEEK_1', 'WEEK_2', 'WEEK_3', 'WEEK_4', 'WEEK_5', 'WEEK_6', 'HOUR_0', 'HOUR_1', 'HOUR_2', 'HOUR_3', 'HOUR_4', 'HOUR_5', 'HOUR_6', 'HOUR_7', 'HOUR_8', 'HOUR_9', 'HOUR_10', 'HOUR_11', 'HOUR_12', 'HOUR_13', 'HOUR_14', 'HOUR_15', 'HOUR_16', 'HOUR_17', 'HOUR_18', 'HOUR_19', 'HOUR_20', 'HOUR_21', 'HOUR_22', 'HOUR_23']

# WebAPIの準備（Tokenの値は上記Webサイトにて取得したものを入力）
headers = {'Authorization': 'Token **************************'}
url = "https://iot.blueomega.jp/spnext/api/v1/weather/forecast.csv"

target_date = ""

# 気象予報値の取得
df_pred = pd.DataFrame()

for city in ["Hiroshima", "Matsue"]:
    
    payload = {
        "point": city
    }
    if len(target_date) > 0:
        payload["date"] = target_date
        
    rs = requests.get(url, headers=headers, params=payload)
    
    df = pd.read_csv(io.StringIO(rs.content.decode('utf-8')), header=0)
    
    # 不具合対応(本来は不要なはず...) BEGIN
    df["flg"] = True
    for i in range(len(df) - 1):
        if df.iloc[i]["forecast_time"] == df.iloc[i+1]["forecast_time"]:
            df.loc[i+1, "flg"] = False
    df = df[df.flg == True]
    # 不具合対応(本来は不要なはず...) END
    
    df.index = pd.to_datetime(df.forecast_time)
    df_pred["TMP_{}".format(city.lower())] = df["temperature"]
    
df_pred["MONTH"] = df_pred.index.month
df_pred["DAY"] = df_pred.index.day
df_pred["WEEK"] = df_pred.index.weekday
df_pred["HOUR"] = df_pred.index.hour

# one-hotエンコーディング
cols = ["MONTH","WEEK","HOUR"]
for col in cols:
    df_pred = df_pred.join(pd.get_dummies(df_pred[col], prefix=col))

# 不足する列を取得
for col in x_cols:
    if col not in df_pred.columns:
        df_pred[col] = 0.0

# 説明変数の取得
X_pred = df_pred[x_cols]
X_pred["rad"] = make_radian_row(pca.transform(X_pred))
X_pred = scaler.transform(X_pred)

df_pred["MW"] = model.predict(X_pred)
```

予測した結果を可視化。
![Unknown.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/114240/af308002-4fa9-58c1-f9ea-f52ad80183d8.png)

1月10日現在、気象予報値を使って予測した12日の電力使用量予測値です。

できた！

12日の電気使用量って本当にこのとおりになるのかな〜。汗