日本最北端の気象予測をしてみた ～来年の楽しい旅行の為に～

はじめに

プログラミングほぼほぼ未経験から、データサイエンティストを目指して学習を開始しました。

まだまだ実践不足ですが、経験を積んで独り立ちを目指します！！

今回の課題を選んだ理由

今年の北海道旅行の模様

毎年、北海道旅行に行っていますが、今年は雪で思うような計画が出来ませんでした。

過去のデータから、翌年の降雪量を予測できれば
対策が講じれるのではないかとの考えから、対応してみました。

実行環境

【OS】
　windows10 64bit
【開発環境】
　Google Colaboratory(Python3.7)
【使用したモジュール】
　・datetime
　・warnings
　・itertools
　・pandas
　・numpy
　・statsmodels
　・matplotlib
　・japanize_matplotlib
　・mean_absolute_error
　・mean_squared_error

アジェンダ

1.　データの収集
2.　japanize-matplotlibのインストール
3.　データの読み込み
4.　データフレームをデータ分析に適した型へ加工
5.　SARIMAモデルを用いた時系列解析
6.　SARIMAモデル検証
7.　予測を実行し、グラフ、数値で可視化
8.　分析結果

1.データの収集

今回は気象庁の気象データを利用させていただきました

国土交通省　気象庁HP

【ダウンロードデータ】
　wakkanai_weather.csv

2.japanize-matplotlibのインストール

python

!pip install japanize-matplotlib

今回は日本語表記が発生する為、japanize-matplotlibをインストール
※□□□□□にならない為の対策です

3.データの読み込み

python

import datetime
import warnings
import itertools
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import japanize_matplotlib
%matplotlib inline


# データの読み込み
# 必要なcoiumnsのみ抽出
weather_df = pd.read_csv("./wakkanai_weather.csv",  engine='python', header=2,encoding = 'shift_jis')
weather_df_temp = weather_df[['年月',  '平均気温(℃)']][2:]
weather_df_snow = weather_df[['年月',  '降雪量合計(cm)']][2:]

ダウンロードした気象データをインポートしました。

取得できるデータは全部取得した為、不要部分は削除

＜今回必要なデータ＞、
・平均気温(℃)
・降雪量合計(cm)
のみ抽出

4.データフレームをデータ分析に適した型へ加工

python

# データの整理
# indexに期間を代入。期間は"1953-01-01"から"2021-12-31"
index = pd.date_range("1953-01-01", "2021-12-31", freq="M")
# "weather_df_temp","weather_df_snow"にindexを代入し、不要なcolumnを削除
weather_df_temp.index = index
weather_df_snow.index = index
del weather_df_temp["年月"]
del weather_df_snow["年月"]

# 検証データとテストデータに分割します。
#  検証データ：「1953/01/01～2009/12/31」
#  テストデータ：「2010/01/01～2021/12/31」
weather_df_temp_tran =weather_df_temp[:684] 
weather_df_snow_tran =weather_df_snow[:684] 
weather_df_temp_test =weather_df_temp[684:] 
weather_df_snow_test =weather_df_snow[684:]

期間をINDEXに代入し、代わりにもともと存在したcolumnを削除（不要の為）
SARIMAモデルの検証用にテストデータを作成

5.SARIMAモデルを用いた時系列解析

今回の解析には
・平均気温
・降雪量
が含まれるものになりますので、季節変動を考慮したSARIMAモデルを適用したいと思います。

python

# orderの最適化関数
def selectparameter(DATA, s):
    p = d = q = range(0, 2)
    pdq = list(itertools.product(p, d, q))
    seasonal_pdq = [(x[0], x[1], x[2], s) for x in list(itertools.product(p, d, q))]
    parameters = []
    BICs = np.array([])
    for param in pdq:
        for param_seasonal in seasonal_pdq:
            try:
                mod = sm.tsa.statespace.SARIMAX(DATA,
                                                order=param,
                                                seasonal_order=param_seasonal)
                results = mod.fit()
                parameters.append([param, param_seasonal, results.bic])
                BICs = np.append(BICs, results.bic)
            except:
                continue
    return parameters[np.argmin(BICs)]

# SARIMAモデルを用いて時系列解析
# 周期は月ごとのデータであることも考慮しs=12
# orderはselectparameter関数の0インデックス, seasonal_orderは1インデックスに格納
#平均気温
best_params = selectparameter(weather_df_temp, 12)
SARIMA_weather_df_temp = sm.tsa.statespace.SARIMAX(weather_df_temp[['平均気温(℃)']], order=best_params[0],
                                             seasonal_order=best_params[1],
                                             enforce_stationarity=False, enforce_invertibility=False).fit()
#降雪量
best_params = selectparameter(weather_df_snow, 12)
SARIMA_weather_df_snow = sm.tsa.statespace.SARIMAX(weather_df_snow[['降雪量合計(cm)']], order=best_params[0],
                                             seasonal_order=best_params[1],
                                             enforce_stationarity=False, enforce_invertibility=False).fit()

SARIMAモデルとは

引用
定常時系列の解析に使われるARMAモデル・SARIMAモデルとは？

6.SARIMAモデル検証

SARIMAモデルの精度を確認していきます。
今回は
・平均絶対誤差(MAE)
・平均平方二乗誤差(RMSE)
を利用して精度評価を行います。

python

#SARIMAモデルの検証

# テストデータと同じ期間で予測
pred_temp_test = SARIMA_weather_df_temp_tran.predict("2010-1-31", "2021-12-31")
pred_snow_test = SARIMA_weather_df_snow_tran.predict("2010-1-31", "2021-12-31")

#今回は
#・平均絶対誤差(MAE)
#・平均平方二乗誤差(RMSE)
#を利用して検証します

from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error

#平均気温
print("平均気温")
#平均絶対誤差(MAE)
print(f'平均絶対誤差(MAE)：{mean_absolute_error(pred_temp_test.values, weather_df_temp_test.values)}')
#平均平方二乗誤差(RMSE)
print(f'平均平方二乗誤差(RMSE):{np.sqrt(mean_squared_error(pred_temp_test.values, weather_df_temp_test.values))}')
print()

#平均気温
print("降雪量合計")
#平均絶対誤差(MAE)
print(f'平均絶対誤差(MAE)：{mean_absolute_error(pred_snow_test.values, weather_df_snow_test.values)}')
#平均平方二乗誤差(RMSE)
print(f'平均平方二乗誤差(RMSE):{np.sqrt(mean_squared_error(pred_snow_test.values, weather_df_snow_test.values))}')
print()

#可視化
fig = plt.figure(figsize=(35, 15))

# グラフ可視化（表示は2019年10月以降）
#予測値は赤色でプロット

#平均気温
ax1 = fig.add_subplot(2, 1, 1)
ax1.set_xlim(datetime.datetime(2010,1,31), datetime.datetime(2021,12,31))
plt.title("平均気温(℃)")
plt.xlabel("年月")
plt.ylabel("平均気温(℃)")
plt.grid(True)
plt.plot(weather_df_temp_test.index, weather_df_temp_test['平均気温(℃)'], marker="o")
plt.plot(pred_temp_test.index, pred_temp_test, color="r", marker="s")

#降雪量
ax2 = fig.add_subplot(2, 1, 2)
ax2.set_xlim(datetime.datetime(2010,1,31), datetime.datetime(2021,12,31))
plt.title("降雪量合計(cm)")
plt.xlabel("年月")
plt.ylabel("降雪量合計(cm)")
plt.grid(True)
plt.plot(weather_df_snow_test.index, weather_df_snow_test['降雪量合計(cm)'], marker="o")
plt.plot(pred_snow_test.index, pred_snow_test, color="r", marker="s")

plt.show()

実行結果

平均気温
平均絶対誤差(MAE)：1.0507855252555975
平均平方二乗誤差(RMSE):1.3376609910211703

降雪量合計
平均絶対誤差(MAE)：14.736495978767644
平均平方二乗誤差(RMSE):26.821912094865493

可視化グラフ

可視化データからわかること
・平均気温は予測値と正解値の誤差はすべての年度において1℃ぐらい
・降雪量合計は予測値と正解値の誤差が大きいことから、「MAE」「RMSE」ともに数値が大きくなっている

上記の理由から、SARIMAモデルは精度は問題ないと推測されます。

7.予測を実行し、グラフ、数値で可視化

それでは予測を開始します。

最初は全期間プロットしました。

すると・・・

になってしまい、細かすぎて全く機能しないグラフが作成されてしましました。


その為、直近3年（2019年10月以降）でプロットするために、

ax.set_xlim(datetime.datetime(2019,10,31), datetime.datetime(2023,4,30))

を追加して、対応しました。

python

# 予測
# 予測期間での予測値を代入
pred_temp = SARIMA_weather_df_temp.predict("2022-1-31", "2023-4-30")
pred_snow = SARIMA_weather_df_snow.predict("2022-1-31", "2023-4-30")

fig = plt.figure(figsize=(35, 15))

# グラフ可視化（表示は2019年10月以降）
#予測値は赤色でプロット

#平均気温
ax1 = fig.add_subplot(2, 1, 1)
ax1.set_xlim(datetime.datetime(2019,10,31), datetime.datetime(2023,4,30))
plt.title("平均気温(℃)")
plt.xlabel("年月")
plt.ylabel("平均気温(℃)")
plt.grid(True)
plt.plot(weather_df_temp.index, weather_df_temp['平均気温(℃)'], marker="o")
plt.plot(pred_temp.index, pred_temp, color="r", marker="s")

#降雪量
ax2 = fig.add_subplot(2, 1, 2)
ax2.set_xlim(datetime.datetime(2019,10,31), datetime.datetime(2023,4,30))
plt.title("降雪量合計(cm)")
plt.xlabel("年月")
plt.ylabel("降雪量合計(cm)")
plt.grid(True)
plt.bar(weather_df_snow.index, weather_df_snow['降雪量合計(cm)'],width=20)
plt.bar(pred_snow.index, pred_snow, color="r",width=20)

plt.show()


#数値データを可視化

#検証データを結合
df_kako = pd.concat([weather_df_temp, weather_df_snow], axis=1)

#予測値を結合
df_pred = pd.concat([pred_temp, pred_snow], axis=1)
df_pred.set_axis(['平均気温(℃)', '降雪量合計(cm)'], 
                    axis='columns', inplace=True)

#検証データと予測値を結合
pd.options.display.precision = 1
df = pd.concat([df_kako, df_pred], axis=0)

print("2019-20年気象データ")
print(df['2019-11-30':'2020-03-31'])
print("\n")
print("2020-21年気象データ")
print(df['2020-11-30':'2021-03-31'])
print("\n")
print("2021-22年気象データ")
print(df['2021-11-30':'2022-03-31'])
print("\n")
print("2022-23年気象データ")
print(df['2022-11-30':'2023-03-31'])

＜平均気温＞

＜降雪量合計＞

＜数値データ＞
※上記グラフを数値化したものになります。

8.分析結果

＜平均気温＞
例年と同じで
・1月～3月にかけて、氷点下を記録
　※本年は4月まで氷点下を記録。予測としては来年も同様な結果となっております。
・最低平均気温はー5度前後を記録
していることが判明しました

＜降雪量合計＞
直近2年では150cm以上の積雪記録をした月が1月は存在しましたが、予測（来年）は存在しません。
但し、2023/01～03にかけては平均的に同じぐらいの積雪が見込まれることが判明しました。

まとめ、及び今後の目標

今回の講座を通じて、データサイエンティストの第一歩を踏み出したが経験不足は否めません。
知識向上には実践あるのみですので
・100本ノックでの実習
・プロボノへの登録
等を積極的に利用して、経験を積んでいきたいと考えています。

目標として、
3年以内に自然言語系に特化したデータサイエンティストとして活躍できるよう、精進します。