過去気象予報を使って太陽光発電電力量を予測する

一般的に日射量と太陽光発電電力量には相関関係があると言われています。
そこで、KCCS APIサービス（以降、KCCS-APIと記載）の過去気象予報を使って太陽光発電電力量を機械学習で予測してみたいと思います。
KCCS-APIの過去気象予報には、予測日射量も含まれます。

|区分|データ種別|内容|
|---|---|---|
|学習データ|説明変数|太陽光発電設備のある緯度経度付近の予測日射量(5月中旬～6月末)|
|学習データ|目的変数|太陽光発電電力量実績(5月中旬～6月末)|
|テストデータ|説明変数|太陽光発電設備のある緯度経度付近の予測日射量(7/1～7/3)|
|テストデータ|目的変数|太陽光発電電力量実績(7/1～7/3)|

なぜ、学習データに、実績日射量ではなく予測日射量を利用するのか、というと、
実際の予測業務では、予測日射量を説明変数とするためです。
明日の太陽光発電量を予測する時点で、明日の実績日射量はないということです。

###過去気象予報を取得してみよう
KCCS-APIでは気象庁の「メソ数値予報」を2016年6月から蓄積しており、ユーザーにて指定いただいた緯度経度をもとに約5kmメッシュ単位で過去の気象予報を一括で取得することができます。
メソ数値予報詳細　https://www.jma.go.jp/jma/kishou/know/whitep/1-3-4.html

例えば、弊社三田事業所の過去気象予報を取得したい場合は、下記のようなCSVファイルを作成し、KCCS-APIの管理画面にジョブを登録するだけで簡単に過去気象予報を取得出来ます。

~~~
KccsMita,35.6916666666667,139.75
~~~


詳細はKCCS-APIサービスHPを確認下さい。
https://www.energy-cloud.jp/specification/#past-weather-forecast

###予測日射量及び太陽光発電電力実績データの取得
弊社社員の住宅に設置された太陽光発電電力量を予測してみることにします。

####予測日射量(過去気象予報データ)の取得
弊社社員の住宅を緯度経度を使って、KCCS-APIから過去気象予報を取得します。
日時フィールドが予測対象の時間となります。1時間ごとにデータがあります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1194651/a0faa822-bc67-1929-4f77-8d9f0b32542f.png)

####太陽光発電電力量実績データの取得
太陽光発電電力量のデータを用意します。
5分毎にデータがあります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1194651/7edb75c8-98f8-5f44-a668-6859afebb39d.png)


###データを整形する
それでは、データを整形していきます。
ここでは、JupyterNotebookを使ってデータ整形、及び分析をしていきたいと思います。

まずはデータを読み込みます。

```py
import pandas as pd

#太陽光発電電力量データの読み込み
df_pv = pd.read_csv("pv.csv")
#過去気象予報データの読み込み
df_forecast = pd.read_csv("WEATHER_FORCAST_RESULT.csv")
df_forecast = df_forecast.sort_values(["日時","予報演算実施日時"])
#必要な情報のみに絞る
df_forecast_ = df_forecast[["日時", "日射量（W/m^2）"]].copy()
df_forecast_ = df_forecast_.rename(columns={'日時':'timestamp', '日射量（W/m^2）':'Insolation'})
df_forecast_['timestamp'] = df_forecast_['timestamp'].astype(str)
~~~

次に時系列を合わせていきます。

```py
#発電電力量は5分ごとのデータのため、1時間ごとのデータに整形
df_pv['timestamp'] = pd.to_datetime(df_pv['timestamp'])
df_pv['timestamp'] = df_pv['timestamp'].dt.round("H")
df_pv = df_pv.groupby('timestamp').sum().reset_index()
~~~

最後に、対象日時でデータをマージします。

```py
#太陽高発電電力量データと過去気象予報データをマージ
df_pv['timestamp'] = df_pv['timestamp'].astype(str)
df_forecast_['timestamp'] = df_forecast_['timestamp'].astype(str)

df_total = pd.merge(df_pv, df_forecast_, on='timestamp', how='inner')
~~~

###太陽光発電電力量を予測する
単回帰分析を使って太陽光発電電力量予測モデルを作成します。

```py
from sklearn.linear_model import LinearRegression

#7/1までのデータを予測モデル作成のデータとして取得します
df_total_ = df_total[df_total["timestamp"] < "2021-07-01" ].reset_index()

x = df_total_[['Insolation']]
y = df_total_[['pv']]

model_lr = LinearRegression()
model_lr.fit(x, y)

plt.plot(x, y, 'o')
plt.plot(x, model_lr.predict(x), linestyle="solid")
plt.xlabel("Insolation")
plt.ylabel("pv")
plt.show()

print('モデル関数の回帰変数 w1: %.3f' %model_lr.coef_)
print('モデル関数の切片 w2: %.3f' %model_lr.intercept_)
print('y= %.3fx + %.3f' % (model_lr.coef_ , model_lr.intercept_))
print('決定係数 R^2： ', model_lr.score(x, y))
~~~

出力結果です。
日射量と発電量に相関関係があることが分かります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1194651/504f729f-29b2-b202-5f4b-79677b019460.png)

作成したモデルを使って、7/1～7/3の予測をしてみます。
青が実績太陽光発電量、オレンジが予測太陽光発電量となります。
細かい誤差はありますが、ある程度は捉えられていることが分かります。

```py
import numpy as np
import matplotlib.pyplot as plt

fig, ax = plt.subplots()

df_total_somedays = df_total[(df_total["timestamp"] >= "2021-07-01") & (df_total["timestamp"] <"2021-07-04")].reset_index()

x_somedays = df_total_somedays[['Insolation']]
y_somedays = df_total_somedays[['pv']]
timestamp_somedays = df_total_somedays[['timestamp']]

#予測を実施
predict_y_somedays = model_lr.predict(x_somedays)

ax.plot(y_somedays, label="measured")
ax.plot(predict_y_somedays, label="predicted")

ax.set_xlabel("time")
ax.set_ylabel("pv")

plt.show()
~~~

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1194651/eb21d603-2dee-ae52-a384-4b380d7bc4c5.png)

###KCCS APIの利用料(ポイント)について

KCCS-APIでの金額の考え方はポイント制となります。
KCCS-APIでは緯度・経度の組み合わせ100地点までは、1日分で100ポイントです。
例えば、1か月分（30日分）のデータを取得すると、1地点でも100地点でも3000ポイントとなります。

安価にデータをダウンロード可能です。
詳細はこちらにてご確認下さい。
https://www.energy-cloud.jp/documents/#service-price-list