2018年7月22日,東京の最高気温は36.5度となり,今年最高を更新した.あまりに暑くてムシャクシャしたので,気象庁から約140年分(1872年1月1日-2018年7月21日)の東京の気温データを入手し,Prophetで分析した.平均・最高・最低気温の全てに関して,1920年付近から上昇し続けており,そのトレンドを考慮してもなお,ここ数日は特に暑いことを確認できた.

なお,分析に用いたNotebookはこちら

環境

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

Prophet

特徴

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

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

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

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

インストール

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

準備

データの入手

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

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

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

area

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

item

最後に,以下の画面で期間を選択する.1872年1月1日から指定可能だったので,選択する.一度にダウンロードできるデータ量に制限があるため,10年ごとに刻んでダウンロードする.

period

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

モジュールのインポート

データの前処理

以下で,詳細を補足する.

codecs.open

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

不要な行・列を除外

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

data

そこで,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日)

分析

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

平均気温

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

scatter_ave

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

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

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

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

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

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

forecast_ave

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

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

components_ave

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

最高気温

同様の分析を行う.詳細は割愛する.

forecast_max

components_min

平均気温と同様の性質が見て取れる.

最低気温

同様の分析を行う.詳細は割愛する.

forecast_min

components_min

平均気温と同様の性質が見て取れる.

感想

地学の知識が中学理科で止まっている私でも,比較的簡単に地球温暖化の様子を確認することができた3.Prophetすごい.

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

  1. デフォルト設定. 

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

  3. あまりに明確にトレンドが出過ぎてちょっと引いた. 

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