schedule2018-07-31

ヒートマップとグラフから札幌の積雪量を分析してみる

前回はProphetを使って札幌の降雪量の分析と予測を行なった。

せっかくなので、別の観点の分析も行う。

ヒートマップを作成

東京の夏が｢昔より断然暑い｣決定的な裏づけのヒートマップが感覚的にわかりやすい。これをめざして見る。

heatmapのカラーはここが参考になる。

# heatmap用DataFrameを生成
years = range(1961, 2018, 1)

depth_each_data = pd.DataFrame()

# シーズン毎に同軸で比較するため
date_1 = datetime.datetime.strptime('2000-11-01', "%Y-%m-%d")

for year in tqdm(years):
    df = all_data.ix[f'{year}-11-01':f'{year+1}-04-30', ['depth']]
    
    # カラム名を年度とする
    df.columns = [f'{year}']
    
    # '2000-11-01'から'2001-04-30' の範囲
    # TODO　閏年の考慮なし。
    arr = range(len(df))
    dｆ.index = [date_1 + datetime.timedelta(days=i) for i in arr]
    
    # date形式に変換
    dｆ.index = df.index.strftime('%y-%m-%d')
    
    # all_dataを更新
    depth_each_data = pd.concat([depth_each_data, df], axis=1)

from matplotlib import dates as mdates
from matplotlib import ticker

fig, ax = plt.subplots(figsize=(18,10)) 
# 1週間ごと
ax.xaxis.set_major_locator(ticker.MultipleLocator(5)) 
sns.heatmap(depth_each_data, cmap="YlGnBu",  vmin=0, vmax=100) # YlGnBu

11月01日から04月30日までの各年度の最深積雪量ヒートマップ。単位はcm。

png

最深積雪量が1００ｃｍを超える年もあるが、見やすさのため上限を1００ｃｍとした。

これを見ると、10年置きに積雪量が多いシーズンを繰り返している。やはり、積雪量には周期性があるようだ。

積雪量は1988年度が最も少なく、1995年度が最も多いことがわかる。

11月にパラパラと降り、根雪（そのシーズン残る雪）は12月中旬くらいから残る。３月末に溶けてからは、積雪はほぼ無い。積雪が多い年は４月上旬まで雪が残っている。

課題は、閏年が考慮されていない点。

# heatmap用DataFrameを生成
years = range(1960, 2018, 1)

amount_each_data = pd.DataFrame()

# シーズン毎に同軸で比較するため
date_1 = datetime.datetime.strptime('2000-11-01', "%Y-%m-%d")

for year in tqdm(years):
    df = all_data.ix[f'{year}-11-01':f'{year+1}-04-20', ['amount']]
    
    # カラム名を変更
    df.columns = [f'{year}']
    
    # '2000-11-01'から'2001-04-30' の範囲
    # TODO　閏年の考慮なし。
    arr = range(len(df))
    dｆ.index = [date_1 + datetime.timedelta(days=i) for i in arr]
    
    # date形式に変換
    dｆ.index = df.index.strftime('%Y-%m-%d')
    
    # all_dataを更新
    amount_each_data = pd.concat([each_data, df], axis=1)

fig, ax = plt.subplots(figsize=(18,10)) 
sns.heatmap(amount_each_data, cmap="YlGnBu",  vmin=0, vmax=30)

11月01日から04月30日までの各年度の降雪量ヒートマップ。単位はcm。

png

11月から4月まで、満遍なく降っている。

降雪量からは、読み取るのが難しい。

シーズン毎の積雪期間と降雪日数をグラフ化して見る

ヒートマップではシーズンの傾向が掴めた。もう少し具体的な変化の様子を見ていく。

from matplotlib import ticker
# 軸の目盛り位置
ax.xaxis.set_major_locator(ticker.MultipleLocator(5))

matplotlibで忘れがちな文法まとめ

# シーズン積雪日数の推移
fig, ax = plt.subplots(figsize=(18,8)) 
plt.rcParams["font.size"] = 18
ax.xaxis.set_major_locator(ticker.MultipleLocator(5)) 
plt.plot(depth_each_data.columns, (depth_each_data > 0).sum(), label="snow_days")

mean = (depth_each_data > 0).sum().mean()
plt.plot(depth_each_data.columns, [mean for _ in range(len(depth_each_data.columns))], label="mean")
plt.plot(depth_each_data.columns, (depth_each_data > 0).sum().rolling(window=5, min_periods=3, center=True).mean(), label="shift_mean")
plt.legend()
plt.show()

png

シーズン毎の積雪期間の日数（snow_days）をグラフ化した。また、平均（mean）とウィンドウサイズが５の移動平均(shift_mean)も計算した。

積雪期間は平均127日。札幌は1年の1/3は雪に覆われていることになる。

移動平均では1978年あたりから約16年周期で変動している。 2009年あたりから積雪が長い期間が続くためであるため、そろそろ雪が少ないシーズンが続くと予想される。

最小２乗法で増減を確認する

from sklearn.linear_model import LinearRegression

x = pd.DataFrame(list(range(len(depth_each_data.columns))))

# 回帰直線モデルを生成
model_lr = LinearRegression()
# 最小二乗法に基づいて、モデル関数を作成
model_lr.fit(x, (depth_each_data > 0).sum())

# シーズン積雪日数の推移
fig, ax = plt.subplots(figsize=(18,8)) 
plt.rcParams["font.size"] = 18
ax.xaxis.set_major_locator(ticker.MultipleLocator(5)) 
plt.plot(x, (depth_each_data > 0).sum(), label="snow_days")
plt.plot(model_lr.predict(x), label="least squares method")

mean = (depth_each_data > 0).sum().mean()
plt.plot(x, [mean for _ in range(len(depth_each_data.columns))], label="mean")
plt.legend()
plt.show()

print('モデル関数の回帰変数 : %.3f' %model_lr.coef_)
print('モデル関数の切片 : %.3f' %model_lr.intercept_)

png

モデル関数の回帰変数 : 0.043
モデル関数の切片 : 125.744

最小２乗法で傾きを求めた。回帰変数0.043と微増する傾向にあった。

降雪日数の変化をグラフ化する

from matplotlib import ticker

# シーズン降雪量のグラフ
fig, ax = plt.subplots(figsize=(18,8)) 
plt.rcParams["font.size"] = 18
ax.xaxis.set_major_locator(ticker.MultipleLocator(5)) 

plt.plot(amount_each_data.columns, (amount_each_data > 0).sum(), label="snowfall_day")
plt.plot(amount_each_data.columns, (amount_each_data >= 5).sum(), label=">= 5cm")
plt.plot(amount_each_data.columns, (amount_each_data >= 10).sum(), label=">= 10cm")
plt.plot(amount_each_data.columns, (amount_each_data >= 20).sum(), label=">= 20cm")
plt.legend()
plt.show()