LoginSignup
250

More than 5 years have passed since last update.

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

Last updated at Posted at 2018-07-23

はじめに

統計開始以来最も暑い夏を生きる皆様,お疲れ様です.あまりに暑くてムシャクシャしたので,気象庁から約140年分(1872年1月1日-2018年7月21日)の東京の気温データを入手し,Prophetで分析しました.

以下は,過去365日分の最高気温の実測値と,2018年7月22日から1ヶ月先までの予測値を表したものです.

2018-07-22-forecast_max.png

また,以下は,最高気温のトレンドと年単位の周期性を表したものです.

2018-07-22-components_max.png

分析の結果,平均・最高・最低気温の全てに関して,1920年付近から上昇し続けており,そのトレンドを考慮してもなお,ここ数日は特に暑いことを確認しました.分析に用いたNotebookはこちらです.

:本記事は,2018年7月22日に個人サイトに投稿した記事を,Qiita向けに再構成したものです.

環境

  • macOS Sierra, 10.12.6
  • Python, 3.6.6
  • Prophet, 0.3

Prophet

特徴

Prophetは,Facebookが開発したオープンソースの時系列解析ライブラリです.Prophetで想定するモデルは:

y(t) = g(t) + s(t) + h(t) + \epsilon(t)

ここで,$y(t)$は時刻$t$における予測値,$g(t)$は時点$t$におけるトレンド成分,$s(t)$は時点$t$における周期成分,$h(t)$は時点$t$における週末などのイレギュラーな成分,そして$\epsilon(t)$は時点$t$における誤差を表します.各成分のパラメータをStanで推定し,モデルを学習します.詳細なモデルは,原論文(Sean J. Taylor and Benjamin Letham, Forecasting at Scale)を参照してください.

Prophetに関する日本語の記事としては,以下のようなものがあります.とても参考になりました.

Prophetの大きな特徴の一つは,時系列解析の前提知識がなくても,誰でも簡単に一定水準の分析をできることです.

インストール

今回はPythonからProphetを使います.pipを使う場合は,pip install fbprophetで,condaを使う場合は,conda install -c conda-forge fbprophetでインストールできます.

準備

データの入手

過去の気象データ・ダウンロード - 気象庁からcsvデータをダウンロードします

  • 地点:東京
  • 項目:
    • 日別値
    • 気温(平均気温,最高気温,最低気温)
  • 期間:1872年1月1日 - 2018年7月21日(53528日)

まず,以下の画面で都道府県および地区を選択します.

2018-07-22-area.png

次に,以下の画面で項目を選択します.降水量や湿度も興味深いですが,一度にダウンロードできるデータ量に制限があるため,今回は気温(平均気温,最高気温,最低気温)のみを選択します.

2018-07-22-item.png

最後に,以下の画面で期間を選択します.一度にダウンロードできるデータ量に制限があるため,10年ごとに刻んでダウンロードします.

2018-07-22-period.png

Downloadフォルダのdata.csvを年代ごとに改名し,data/raw/以下に次のように保存します.

  • data_1870.csv
  • data_1880.csv
  • data_1890.csv
  • data_1900.csv
  • data_1910.csv
  • data_1920.csv
  • data_1930.csv
  • data_1940.csv
  • data_1950.csv
  • data_1960.csv
  • data_1970.csv
  • data_1980.csv
  • data_1990.csv
  • data_2000.csv
  • data_2010.csv

モジュールのインポート

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

sns.set_style(style='ticks') # グラフスタイルの指定

データの前処理

years = range(1870, 2020, 10)
read_path = 'data/raw/'

all_data = pd.DataFrame()

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

    # 補足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, 4, 7])

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

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

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

以下で,詳細を補足します.

codecs.open

Shift-JISで作成されているため,普通にpandas.read_csv()を使っても,ファイルを読み込めません.pandasでread_csv時にUnicodeDecodeErrorが起きた時の対処 (pd.read_table()) - Qiitaを参考に,codecs.open()を利用します.

2018年7月25日追記:@skotaroさんのご指摘の通り,pandas.read_csv(encoding='Shift-JIS')で普通に読み込めました.私の不勉強でした.ご指摘ありがとうございました!

不要な行・列を除外

元データをそのまま読み込むと,以下のように不要な行・列が含まれてしまいます.

2018-07-22-data.png

そこで,pandas.read_table()のパラメータskiprowsusecolsを用いて,分析に必要な部分のみ抽出します.

すべてnp.nanの行を削除

気象庁のサイト上では1872年1月1日から指定可能でしたが,最初の数年分はデータが含まれていませんでした.pandas.DataFrame.dropna(how='all')で,すべての値がnp.nanの行を削除します.これにより全体の約2%にあたる1272日分のデータが除外されます.

  • 地点:東京
  • 項目:
    • 日別値
    • 気温(平均気温,最高気温,最低気温)
  • 期間:1875年6月10日 - 2018年7月21日(52266日)

分析

平均気温,最高気温,および最低気温の分析を行します.

平均気温

まず,すべての点を散布図でプロットしてみます.

plt.figure(figsize=(8, 6))
plt.plot(all_data.index, all_data.average, 'o', alpha=.1)
plt.xlabel('Date')
plt.ylabel('Average temparature')
plt.show()

2018-07-22-scatter_ave.png

140年間で平均気温がじわじわ上がっていく様子がわかります.Prophetでより詳細に分析します.

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

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

# 学習
m.fit(ave_data)

Prophetで時系列解析を行うためには,以下のカラムを持つpandas.DataFrameを生成する必要があります.

  • ds:日付カラム.
  • y:予測対象カラム.

また,Prophetではデフォルトでyearly_seasonalityweekly_seasonalityを仮定して学習します.気温は曜日に依存しないと考えられるため,weekly_seasonalityを除外します.

以上で学習したモデルmを使い,将来30日間の平均気温を予測します.以下のように,m.make_future_dataframeDataFrameを生成すると楽です.

# 将来30日間を予測
future = m.make_future_dataframe(periods=30)
forecast = m.predict(future)

予測結果をm.plot()でプロットします.約140年分プロットするとわけがわからないので,直近1年分のみを表示するようにします.

m.plot(forecast)
plt.xlim(future.ds.iloc[-365-30], future.ds.iloc[-1])
plt.xlabel('Date')
plt.ylabel('Average')

2018-07-22-forecast_ave.png

黒点が実測値,青線の実線が予測値,そして薄い青の塗りつぶしが80%信頼区間1を示します.ここ数日の平均気温は,信頼区間の上限ギリギリに張り付いていることが確認できます.例年から統計的に予測したとき,いかにここ数日の気温が高いかわかります2.また,例年の傾向を考慮すると,気温のピークは8月初旬であるという絶望的な事実がわかります.

次に,m.plot_components(forecast)で周期性とトレンドを確認します.

m.plot_components(forecast)

2018-07-22-components_ave.png

上の図がトレンド,下の図が年内の周期性を表します.東京の平均気温は1920年あたりから上がり続けていることがわかります.また,年内で見ると,気温のピークは7月下旬から8月上旬にあることがわかります.

最高気温

同様の分析を行います.詳細は割愛しますが,平均気温と同様の性質が確認できます.

2018-07-22-forecast_max.png

2018-07-22-components_max.png

最低気温

同様の分析を行います.詳細は割愛しますが,平均気温と同様の性質が確認できます.

2018-07-22-forecast_min.png

2018-07-22-components_min.png

おわりに

地学の知識が中学理科で止まっている私でも,今年の夏はヤバイということがわかりました3.無理は禁物です.ご自愛くださいませ.

あと,この記事を書き終わるころ,東洋経済さんの素晴らしい記事を見つけました(東京の夏が「昔より断然暑い」決定的な裏づけ: 過去140年の最高気温をビジュアル化 - 東洋経済オンライン4.こんなにわかりやすいヒートマップを初めて見ました.いつかこんなビジュアライゼーションをしてみたいです.

最後まで読んでくださり,ありがとうございました!


  1. デフォルト設定. 

  2. 信頼区間を外れた異常値は他にも結構存在しますが,連日,しかも同じ方向に外れることは珍しいように思います. 

  3. あまりに明確に温暖化のトレンドが出過ぎてちょっと引きました. 

  4. 一応,1920年ごろから気温が上がり始めるという分析結果は,本記事と整合しています.良かった. 

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
250