schedule2018-07-31

Prophetで札幌の積雪量の推移を分析してどれだけ積もるか予測する

趣旨

  • Prophetがすごい
  • 札幌の降雪量で確認する
  • 別の方法でも検証してみる

はじめに

あまりに暑いので,140年分の気温をProphetで分析した

この記事を参考に、自分も分析してみたくなり始めました。

内地(北海道の人は本州をこう呼ぶ)は暑いようですが、札幌のここ一週間の最高気温は32℃です(yahoo調べ)。夜、窓を開けていると寒いぐらいです。

気温はいろんな方が分析されているので、季節を先取り積雪量について分析してみました。 地球温暖化の影響もあって積雪量は少なくなっていく印象ですが、実際はどうなのでしょうか。 札幌のこれまでの積雪量の推移来シーズン(2018-2019)の積雪量の予測を検証していきます。

Prophetの解析は、@haltaroさんの手順通りです。 詳しい内容は、@haltaroさんの記事をご参照ください。 また、全く同じでは申し訳ないので、Prophetの結果をヒートマップなど別の角度から検証します。

以下は、約60年間のデータから来シーズンの札幌の降雪量をProphetで予測した結果です。

png

以下は、そのトレンドと年単位の周期です。

png

以下は、積雪量のヒートマップです。

png

環境

  • mac OS High Sierra, 10.13.4
  • Python 3.7.0
  • Prophet 0.3

データ

同じく、データは過去の気象データ・ダウンロード - 気象庁から入手しました。

1961年01月01日から最深積雪量降雪量の日合計の観測記録が残っています。

  • 地点:札幌
  • 期間:1961年01月01日 - 2018年07月29日 (21029日)
  • 項目:
    • 最深積雪量(cm)
    • 降雪量の日合計(cm)

参考の通りに、10年毎の記録を入手しました。

モジュールのインポートと前処理

import numpy as np # 数値計算
import pandas as pd # DataFrame
import matplotlib.pyplot as plt # グラフ描画
import seaborn as sns # グラフ描画設定
from tqdm import tqdm # forループの進捗状況確認
import codecs # Shift-JIS読み込み用
from fbprophet import Prophet # prophet

import datetime

sns.set_style(style='ticks') # グラフスタイルの指定
years = range(1960, 2020, 10)
read_path = 'snow_sapporo/'

all_data = pd.DataFrame()

for year in tqdm(years):
    data_file = f'data_{year}.csv'

    # 補足1:普通にpd.read_csvすると読み込めない
    with codecs.open(read_path+data_file, "r", "Shift-JIS", "ignore") as f:
        # 補足2:先頭数行のデータは使用しない.
        # また,不要なカラムは除外して利用する.
        data = pd.read_table(
            f, delimiter=",", skiprows=5, index_col=0, 
            usecols=[0, 1, 5])

    # datetime形式に変換
    data.index = pd.to_datetime(data.index)
    # カラム名を変更
    data.columns = ['depth', 'amount']

    # all_dataを更新
    all_data = pd.concat([all_data, data])

# 補足3:すべてNanの行は削除
all_data = all_data.dropna(how='all')

cファイル名をdata_{year}にしただけで、全て読み込めました。

雪積量の推移

全ての点を散布図でプロットした結果。

plt.figure(figsize=(8, 6))
plt.plot(all_data.index, all_data.depth, 'o', alpha=.1)
plt.xlabel('Date')
plt.ylabel('depth_of_snowfall (cm)')
plt.show()

最深積雪量の散布図

png

数年置きに雪積量が上下している。 全体の傾向としては、積雪が減少している訳ではなさそう。

plt.figure(figsize=(8, 6))
plt.plot(all_data.index, all_data.amount, 'o', alpha=.1)
plt.xlabel('Date')
plt.ylabel('amount_of_snowfall (cm)')
plt.show()

降雪量の日合計の散布図

png

大きな変化は見られない。

2012年辺りの最深積雪量は高いのに、1日の降雪量はそこまで大きくない。気温が低く、溶けずに積もっていったと思われる。

Prophetで来シーズンの積雪を予想する

2018-2019シーズンの積雪量を予想します。

# Prophet用DataFrameを生成
depth_data = pd.DataFrame()
depth_data['ds'] = all_data.index
depth_data['y'] = all_data.depth.values

# Prophetモデルの構築
m = Prophet(weekly_seasonality=False)

# 学習
m.fit(depth_data)
# 将来9ヶ月を予測
future = m.make_future_dataframe(periods=270)
forecast = m.predict(future)
m.plot(forecast)
plt.xlim(future.ds.iloc[-910], future.ds.iloc[-1])
plt.xticks(rotation=90)
plt.xlabel('Date')
plt.ylabel('depth_of_snowfall (cm)')
plt.show()

積雪量のProphetの結果です。

png

黒点が実測値,青線の実線が予測値,そして薄い青の塗りつぶしが80%信頼区間を示します。

今年も12月に積もり始め、例年通りの降雪量になりそうです。

2016年12月は連日「うんざり」 札幌50年ぶり96センチ 12月との報道があり、12月としては1966年以来の大雪だそうです。信頼区間から大きくはみ出していますね。

積雪量は半波整流波のようなモデルですが、上手く予測することができました。

m.plot(forecast)
plt.xlim(future.ds.iloc[-4000], future.ds.iloc[-1])
plt.xticks(rotation=90)
plt.xlabel('Date')
plt.ylabel('depth_of_snowfall (cm)')
plt.show()

png

積雪量がない期間の80%信頼区間が常に広いです。 0に収束すると考えていましたが、違うようですね。 0付近か、変化のない期間は苦手なのでしょうか。

失敗例として、各年度の10月から4月だけのデータから予測した結果も載せておく。 シーズンのデータだけで良いと思ってました。

失敗

データがない5月から9月までの期間の予想が荒ぶっている。 データが揃っている場合は、大きな欠損期間はない方が良さそう。

予測モデルの構成要素

トレンドと周期性についても分析結果を見てみます。

#トレンドや周期性を表示する
m.plot_components(forecast)

png

トレンドはどちらかというと右肩上がり。 積雪量は減少せず、むしろ増えていっている。

15年周期で変動しているようにみえる。ここ10年程は上昇しているため、そろそろ下降に転じそう。

年周期性は11月から積もり始め、2月下旬が最も高く4月頃には雪が溶けている。 緩やかに積もり、急激に溶けていく。

積雪がない時期がマイナスになっているのはあっているのだろうか?

ここまでのまとめ

  • 積雪量は減少していない。
    • むしろ、増加している。
  • 積雪量は約15年の周期で変動している。
    • ここ10年程は増加傾向にあるため、来シーズンから減少に転じると思われる。
  • 来年も例年通り、12月から3月末までの積雪となる。

Prophet以外の方法でも、同様のことが言えるか確かめる。

こちらに続く。