【文系・未経験者】日本がハワイ化するのは何年後か予想してみた

温暖化の影響で日本の気候は熱帯気候に近づいているという。リゾートの代名詞「ハワイ」に平均気温が追いつくのは今から何年後の話しになるのかを予想してみました。

開発の環境

使用言語：Python3
使用機種：MacBook Pro
使用ブラウザ：Google Chrome
開発環境：Google Colaboratory

自身の状況

Aidemy Premiumにて「データ分析6ヶ月コース」を受講。
完全にプログラミング初心者で、プログラムの「プ」の字も理解できていなかったところからスタートしています。

やりたいこと

現時点のハワイの年間平均気温を取得し、その気温をベンチマーク。
日本（東京）の年々上がっている平均気温を過去のデータから今後のデータを推測し、現時点のハワイの平均気温にいつ頃追いつくのかの予測を出したいと思います。

作業1.データの収集

初めに各種必要となるデータを集めていきます。
今回必要なデータはハワイの平均気温データと、日本（東京）の平均気温データです。
データは気象庁のサイト(https://www.jma.go.jp/jma/menu/menureport.html)にて提供があったため、そこから拝借します。

■データ1.東京の平均気温データ
サイトを見ると、1872年頃からのデータが閲覧できました。
「過去の気象データ・ダウンロード」のページからデータをCSV形式でダウンロードできるようなので、地域を「東京」、項目を「月平均気温」に設定し、ダウンロードします。

■データ2.ハワイの平均気温データ
「地点別データ・グラフ (世界の天候データツール（ClimatView 月統計値）)」のページより、世界の天候のデータもダウンロードできました。
地域を「オアフ島」を選択、「グラフ表示期間」をすべての年にしてみます。

作業2.データの確認

ひとまず、ダウンロードしたデータをテキストエディタで見てみます。

data.csv（東京のデータ）

climat_91285_198206-202304.csv（ハワイのデータ）

双方で年月の分け方が違ったり、使わないであろうデータが入っていたりとするので、読み込み後に整形が必要そうなデータです。

作業3.データの読み込み

ダウンロードしてきたデータを、Google Colaboratoryにて読み込みます。
作業.2で確認した際に、日本語での記述があったため、文字コードは"shift_jis"を指定しておきます。また、必要なデータは、年月と平均気温だけなので、csvの読み込み時に、"usecols"を使い、整理してしまいます。

from google.colab import files
uploaded = files.upload()

↑この技術はすごいなぁ、といつも感嘆します。
さて、アップロードしたファイルをちゃんとアップロードできているかheadで確認してみます。

import pandas as pd
tokyo_temp = pd.read_csv('./data.csv',  engine='python', header=1,index_col=0,encoding = 'shift_jis',usecols=[0,1])
hawaii_temp = pd.read_csv('./climat_91285_198206-202304.csv',  engine='python', header=1,index_col=0,encoding = 'shift_jis',usecols=[0,1,2])
print(tokyo_temp.head())
print()
print(hawaii_temp.head())

head結果

ちゃんと読み込みができているようです。
受講したての頃はこれでも嬉しかった記憶がありますが、ここからデータの加工を始めていきます。

作業4.データの加工

東京のデータは年月がひとつにまとまっており、ハワイのデータは別々となっているため、ここから合わせていきます。具体的には、"/"で分割し、不要な行を削除していきます。項目名も「年」「月」「気温」で統一します。多分もっと綺麗なコードがあるんだと思いますが、いろいろ調べていき、ここが限界でした。

base_file = "data.csv"
use_file = "data2.csv"

with open(base_file, "r", encoding="shift_jis") as cr:
  lines = cr.readlines()
  lines = map(lambda v: v.replace('/', ','), lines)
  result = "".join(lines).strip()

with open(use_file, "w", encoding="shift_jis") as cw:
  cw.write(result)

tokyo_temp = pd.read_csv('./data2.csv',  engine='python', header=1, index_col=0, encoding = 'shift_jis', usecols=[0,1,2], names=["年","月","気温"])
tokyo_temp = tokyo_temp.reset_index()
tokyo_temp = tokyo_temp.drop([0,1], axis=0)

print(tokyo_temp.head())
print()
print(hawaii_temp.head())

出力結果

さらに、年平均の気温とするため、pandasのgroupbyで年毎に平均（mean())でまとめていきます。
が、ここで想定外の大苦戦となりました。
1つ目が「年」でまとめたことでグループラベルがインデックスになってしまい、混乱。2つ目がハワイの年平均を出そうとしたところ、うまく計算ができない、と言う問題が発生しました。（※全て23.58…という気温になってしまう）

1つ目はgroupbyの処理の際に、as_index=Falseを指定することで解決。
2つ目は原因が特定できずに時間がかかりましたが、データにテキストの部分があることがわかったため、floatに変換しようとするも失敗。
エラーメッセージをよく読むと、データ末尾にスペースが入っている箇所があり、それが邪魔をしていたことがわかりました。
先ほどの"/"を削除した要領でスペースを削除し、groupbyをかけていきます。

with open(base_file, "r", encoding="shift_jis") as cr:
  lines = cr.readlines()
  lines = ["年,月,気温,品質,均質"] + lines[4:]
  lines = map(lambda v: v.replace('/', ','), lines)
  result = "".join(lines).strip()

with open(use_file, "w", encoding="utf-8") as cw:
  cw.write(result)

# ハワイの処理用に追加
with open(base_file2, "r", encoding="shift_jis") as cr:
  lines = cr.readlines()
  lines = ["年,月,気温,品質,均質"] + lines[4:]
  lines = map(lambda v: v.replace(' ',''), lines)
  result = "".join(lines).strip()

with open(use_file2, "w", encoding="utf-8") as cw:
  cw.write(result)

#groupbyで年平均気温算出
tokyo_temp = pd.read_csv('./data2.csv', encoding = 'utf-8')
tokyo_temp = tokyo_temp.reset_index()
hawaii_temp = pd.read_csv('./climat_91285_198206-202304-2.csv',  engine='python',  encoding='utf-8')
hawaii_temp = hawaii_temp.reset_index()
df = tokyo_temp.groupby(tokyo_temp['年'], as_index=False).mean()
df2 = hawaii_temp.groupby(hawaii_temp['年'], as_index=False).mean()

ここまでを一旦グラフに表示してみたいと思います。
グラフ表示のモジュールをインポート。

import matplotlib.pyplot as plt
%matplotlib inline

data_x = df['年']
data_y = df['気温']

data_x2 = df2['年']
data_y2 = df2['気温']

fig, ax = plt.subplots()

ax.plot(data_x, data_y, linestyle='solid', marker='o')
ax.plot(data_x2, data_y2, linestyle='solid', marker='o')

ax.set_xlabel('year')
ax.set_ylabel('temp')

plt.show()

グラフ化。東京（青）、ハワイ（オレンジ）、ハワイあったかい。

東京の平均気温はここ140年で上昇傾向が見受けられます。ハワイはここ40年でそこまでの変動は見られません。
なので、ハワイはこの40年の平均気温を算出し、東京はあと何年後にその気温にたどり着くのか？を算出していきたいと思います。また、東京の気温のデータ、最初と最後は1年分のデータが入っていないため、削除します。

hawaiimeantemp = df2.mean()
print(hawaiimeantemp['気温'])

>>23.547965367965368

算出できました。これを基準とし、未来の東京の平均気温を時系列分析で予想していきたいと思います。コードもぐちゃぐちゃしてきたので整理します。

#東京の気温データの読み込みと不要箇所の削除
tokyo_temp = pd.read_csv('./data2.csv', encoding='utf-8')
tokyo_temp = tokyo_temp.reset_index()
hawaii_temp = pd.read_csv('./climat_91285_198206-202304-2.csv', engine='python', encoding='utf-8')
hawaii_temp = hawaii_temp.reset_index()
df = tokyo_temp.groupby(tokyo_temp['年'], as_index=False).mean()
df = df.drop(df.iloc[:4].index)#1874年まで1年分のデータがないため削除
df = df.drop(df.iloc[147:].index)#2023年は1年分のデータがないため削除
df['年'] = pd.to_datetime(df['年'], format='%Y')
df.set_index('年', inplace=True)

data_x = df.index
data_y = df['気温']

#ハワイも同様に不要箇所の削除と「ハワイ平均気温」（今回の目標値）の決定
df2 = hawaii_temp.groupby(hawaii_temp['年'], as_index=False).mean()
df2 = df2.drop(df2.iloc[41:].index)  # 2023年は1年分のデータがないため削除
df2['年'] = pd.to_datetime(df2['年'], format='%Y')
df2.set_index('年', inplace=True)
data_y2 = df2['気温']

hawaiimeantemp = df2.mean()
print(hawaiimeantemp['気温'])

作業5.時系列分析（SARIMAモデル）

データの準備が整ったので、いよいよ分析を行っていきます。
今回、未来の数値を予測したいことと、勉強中の添削課題にもなっていたことで多少は見慣れているSARIMAモデルを使ってみます。

import statsmodels.api as sm
start_year = 2022
end_year = 2262  # 補足あり

model = sm.tsa.SARIMAX(data_y, order=(0, 1, 1), seasonal_order=(1, 0, 1, 12))
model_fit = model.fit(disp=0)
pred_japan = model_fit.predict(start=str(start_year), end=str(end_year), dynamic=True)

まず、分析にあたり、「end_year」の設定に四苦八苦しました。
当初は2100年としていましたが、全然目標値に届かない状態だったため、
思い切って3000年まで伸ばしたところエラーが発生。

エラーコード

要は3000年なんてないよ、というエラーでした。じゃあいつまでならできるのよ？と散々調べたところ、どうやらpandasのDatetimeは2262-04-11までしかサポートしていない様子。そのため、2262年を「end_year」に設定しました。

次のハードルはSARIMAモデルのパラメータ決定の部分でした。学習した内容を何度も読み返し、AICの値がもっとも低いものを採用してみることとしました。学習した際はBICの値が低いものを採用する関数でした。

# 各パラメータ範囲
p = d = q = range(0, 3)
sp = sd = sq = range(0, 2)
 
# p, d, q　の組み合わせリストの作成
pdq = list(itertools.product(p, d, q))
 
# P, D, Q　の組み合わせを列挙するリストを作成
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(sp, sd, sq))]
 
best_result = [0, 0, 10000000]
for param in pdq:
    for param_seasonal in seasonal_pdq:
        try:
            mod = sm.tsa.SARIMAX(data_y,
                          order = param,
                          seasonal_order=param_seasonal,
                          enforce_stationarity=True,
                          enforce_invertibility=True)
             results = mod.fit()
            print('order{}, s_order{} - AIC: {}'.format(param, param_seasonal, results.aic))
 
            if results.aic < best_result[2]:
                best_result = [param, param_seasonal, results.aic]
        except:
            continue
             
print('BESTAIC:', best_result)

いろいろな方の解説ブログやドキュメント参考に作っていきます。ここはもう一度深く勉強が必要だと痛感しました。

BEST AIC

パラメータ数値がわかったところで、その数値を入れたコードを書き、モデル化していきます。

model = sm.tsa.SARIMAX(data_y, order=(0, 1, 1), seasonal_order=(1, 0, 1, 12))
model_fit = model.fit()
pred_japan = model_fit.predict(start=str(start_year), end=str(end_year))

モデル化できたら順を追ってグラフにしてみます。
まずは日本の平均気温と、目標とするハワイの平均気温です。
ハワイの平均気温は到達地点とわかるように、点線で水平に引きます。
axhlineを使用すると実現できることがわかり、採用します。

# グラフの表示
plt.plot(data_x, data_y, label='Japan temp')　#過去から今までの日本の気温
plt.axhline(y=hawaiimeantemp['気温'], color='r', linestyle='--', label='hawaii temp')  #目標とするハワイの平均気温
plt.xlabel('年')
plt.ylabel('気温')
plt.legend()
plt.show()

過去から今までの日本の気温（青）と目標とするハワイの気温（赤）

次に現在〜未来の日本の平均気温予測をグラフ化します。

# グラフの表示2
plt.plot(pred_japan.index, pred_japan, label='Japan temp future')#未来の日本の気温予想
plt.xlabel('year')
plt.ylabel('temp')
plt.legend()
plt.show()

これから先の日本の予想年間平均気温

それぞれがグラフ化できたところで、1つのグラフに表示していきます。

# グラフの表示
plt.plot(data_x, data_y, label='Japan temp') #過去から今までの日本の気温
plt.axhline(y=hawaiimeantemp['気温'], color='r', linestyle='--', label='hawaii temp')  #目標とするハワイの平均気温
plt.plot(pred_japan.index, pred_japan, label='Japan temp future')#未来の日本の気温予想
plt.xlabel('year')
plt.ylabel('temp')
plt.legend()
plt.show()

日本のこれからの年間平均気温予測

着実に気温は上昇していくことが予測されますが、このモデルでは2262年でもまだまだハワイの気温には届かない予測となりました。
日本がハワイになるのを待つより、ハワイに行ってしまった方が良さそうです。

まとめ

pandasのDatetimeの制約により、「ハワイに到達するのはいつか？」というのは出すことができませんでした。（おそらく、出す方法はいくらでもあるのだとは思います）
今回、このモデルを組むにあたり、間違いなくいい経験とはなりましたが、SARIMAモデルのパラメータの部分等、学び直しの必要や、もっと深掘りして学ぶ必要があることが改めてわかりました。公式ドキュメントや、他の方の解説ブログ等を読み漁りましたが、まだまだ自分の知識として吸収できているのか？と問われると自信がありません。

今後は、他のモデルでも同様の予測を出してみたり、pandasの制約をクリアする方法も模索してみたいと思います。

お見苦しい点は多々あったかと思いますが、最後まで素人にお付き合いいただき、ありがとうございました。