温室効果ガスによる平均気温の上昇予測

気象庁の温室効果ガス濃度データと平均気温を利用して、今後の気温上昇を予測することを目的とします。

その際、二酸化炭素濃度とメタンガスにおける気象庁データの都合より、与那国島でのデータを用いて予測を行います。

環境

Python : 3.9.7 64bit
MacBook Pro  OS: macOS Monterey
Visual Studio Code  バージョン: 1.62.0 (Universal)

データセット

与那国島　月平均気温 （1987年1月〜2021年10月）
与那国島　月平均二酸化炭素濃度 （1997年1月〜2021年6月）
与那国島　月平均メタンガス濃度 （1998年1月〜2021年6月）

CSVデータの確認と前処理

まずデータをテキストエディットより開き確認します。　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　【月平均気温】　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　年月,平均気温（℃）,品質情報,均質番号

【月平均二酸化炭素濃度　月平均メタンガス濃度】　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　年,月,月平均濃度（ppm）,

月平均気温データでは、年月が一緒になっているので、CSVファイルを成形し直します。またCSVファイルは日本語が混ざっているので、開く際はencoding="Shift_JISで読み込みます。以下、コードと成形後のデータです。

# データの整形
in_file = "yo_tempe_data.csv"
out_file = "yo_tempe_data1.csv"

# CSVファイルを読み込む
with open(in_file, "r", encoding="Shift_JIS") as fr:
   lines = fr.readlines()

# ヘッダーを付け替える
lines = ["年,月,平均気温,品質,均質\n"] + lines[1:]
lines = map(lambda v: v.replace('/', ','), lines)
result = "".join(lines).strip()

# 結果をファイルへ出力
with open(out_file, "w",encoding = 'shift_jis') as fw:
   fw.write(result)

これで、平均気温とその他データとカラムを合わせることができました。

CSVファイルの読み出し

必要な情報は年,月,平均気温（℃）なので、usecols=[0,1,2]で品質情報,均質番号を読み込み時点で省きます。

#与那国島 月平均気温 CSVの読み込み (読み込みと同時に必要な列のみ取り出し)
ave_temp = pd.read_csv("yo_tempe_data1.csv",encoding = 'shift_jis',usecols=[0,1,2])

#与那国島 月平均二酸化炭素濃度 CSVの読み込み (読み込みと同時に必要な列のみ取り出し)
ave_co2 = pd.read_csv("co2_monthave_yon.csv",encoding = 'shift_jis',usecols=[0,1,2])

#与那国島 月平均メタンガス濃度 CSVの読み込み (読み込みと同時に必要な列のみ取り出し)
ave_ch4 = pd.read_csv("ch4_monthave_yon.csv",encoding = 'shift_jis',usecols=[0,1,2])

データの前処理

月平均気温については（1987年1月〜2021年10月）までのデータがありますが、月平均二酸化炭素濃度 については（1997年1月〜2021年6月）、月平均メタンガス濃度 については（1998年1月〜2021年6月）のデータとなりますので、月平均メタンガス濃度のデータに合わせて他データを処理します。

#月平均気温 を1996年12月まで行を削除する（二酸化炭素濃度と条件を合わせるため）
ave_temp = ave_temp.drop(range(0,120))
ave_temp = ave_temp.drop(range(402,418))

#pd.set_option('display.max_rows', 300)
#print(ave_temp)

#二酸化炭素濃度データの末尾に不要な文章があるため行を削除する
ave_co2 = ave_co2.drop(ave_co2.index[[294,295,296,297]])


#メタンガス濃度データの末尾に不要な文章があるため行を削除する
ave_ch4 = ave_ch4.drop(ave_ch4.index[[282,283,284,285]])


#二酸化炭素濃度データとメタンデータの2020年6月以降の行を削除する（メタン濃度のデータがないため）
ave_co2 = ave_co2.drop(range(282,294))
#pd.set_option('display.max_rows', 300)
#print(ave_co2)

欠損値を以下メソッドで確認します。

# 欠損値の数
print(ave_temp.isnull().sum())


# 欠損値の数
print(ave_co2.isnull().sum())


# 欠損値の数
print(ave_ch4.isnull().sum())

年 0
月 0
平均気温 0
dtype: int64
年 0
月 0
二酸化炭素濃度の月平均値(与那国島)[ppm] 0
dtype: int64
年 0
月 0
メタン濃度の月平均値(与那国島)[ppb] 0
dtype: int64

と出力されますが、実際のデータには欠損値があり、『NaN』ではなく『---』となっているため、検出できないようです。

実際にデータを目視し欠損値を確認後、直接欠損値の補完をfillna()メソッドを利用して対応します。

#二酸化炭素濃度データ欠損値の補完
ave_co2.iloc[225, 2] = NA
ave_co2 = ave_co2.fillna(method="ffill")

#pd.set_option('display.max_rows', 300)
#print(ave_co2)


#メタンガス濃度データ欠損値の補完
ave_ch4.iloc[213, 2] = NA
ave_ch4.iloc[12, 2] = NA
ave_ch4 = ave_ch4.fillna(method="ffill")

#pd.set_option('display.max_rows', 300)
#print(ave_ch4)

基本統計量を確認します。

ave_ch4 = ave_ch4.astype(np.float64)　#文字列になっているので変換
ave_co2 = ave_co2.astype(np.float64)　#文字列になっているので変換

# 基本統計量
print(ave_temp.describe(include='all'))
print(ave_co2.describe(include='all'))
print(ave_ch4.describe(include='all'))

以下出力結果です。

後の学習を行いやすくするために、二酸化炭素濃度データにメタンガスデータを追加します。

#CO2データにメタンデータを追加する
ave_co2["メタン濃度の月平均値(与那国島)[ppb]"] = ave_ch4["メタン濃度の月平均値(与那国島)[ppb]"]

データの差分を取る

回帰分析を利用して分析するため、隣り合う2時点間での相関を解消します。今回扱うデータは月平均気温と月平均温室効果ガス濃度です。　　　　　　　　　　月平均気温は季節により変動しますし、月平均温室効果ガス濃度も1年間の間で周期的な変動があるかもしれません。よって、1年間の間で周期変動するデータ（季節変動）があることを想定し、【当月データ】から【前月データ】の値を引く（階差系列）ことでトレンドを取り除いた月平均気温変化分と月平均ガス濃度変化分とします。　　　　　　　　　　　　　　　　　　　　同時に、先頭データは前月が存在せず欠損値となってしまうため、インデックス0の欠損値も処理します。

#各データの差分を計算
diff_temp = ave_temp["平均気温"].diff().fillna(method='bfill') #インデックス0の欠損値についても処理する
#print(diff_temp)


diff_warm_gass = ave_co2[["二酸化炭素濃度の月平均値(与那国島)[ppm]","メタン濃度の月平均値(与那国島)[ppb]"]].diff().fillna(method='bfill') #インデックス0の欠損値についても処理する

学習用データとテスト用データに分ける

目的変数を平均気温(diff_tempに格納)、説明変数を二酸化炭素濃度とメタンガス濃度の2データ(diff_warm_gassに格納)とし、それぞれを2:8でテスト用データと学習用データに分けます。

学習用データは1997年1月〜2015年9月まで期間、テストデータは2015年10月〜2020年6月までの期間となっています。

# データを学習用と予測用に分離する
Y = diff_temp
X = diff_warm_gass

X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2, random_state=0, shuffle=False)
#print(X)

モデルの選択

どのモデルが良いかを線形回帰・ラッソ回帰・リッジ回帰から選択します。

max_score = 0
best_model = ""
# 線形回帰
model = LinearRegression()
model.fit(X_train, Y_train)
score = model.score(X_test, Y_test)
if max_score < score:
   max_score = score
   best_model = "線形回帰"
# ラッソ回帰
model = Lasso()
model.fit(X_train, Y_train)
score = model.score(X_test, Y_test)
if max_score < score:
   max_score = score
   best_model = "ラッソ回帰"
# リッジ回帰
model = Ridge()
model.fit(X_train, Y_train)
score = model.score(X_test, Y_test)
if max_score < score:
   max_score = score
   best_model = "リッジ回帰"

print("モデル:{}".format(best_model))
print("決定係数:{}".format(max_score))

print('決定係数(train):{:.3f}'.format(model.score(X_train,Y_train)))
print('決定係数(test):{:.3f}'.format(model.score(X_test,Y_test)))

以下出力結果です。

モデル:リッジ回帰
決定係数:0.38540234497824344
決定係数(train):0.475
決定係数(test):0.385

よってモデルはリッジ回帰を利用しますが、決定係数が低い値です。

#モデル作成
model = Ridge()
model.fit(X_train, Y_train)

また、この時の傾きと切片を確認します。

coef = list(model.coef_) #傾き
intercept = model.intercept_ #切片

print(coef)
print(intercept)

以下、各値が出力されました。
傾き：[-0.393546658686044, -0.026502783665514264]
切片：0.17738695202156465

予測

現状怪しい学習結果ではありますが、どのような予測になるのか確認します。トレーニングデータとテストデータをそれぞれ予測させ、その結果を結合します。

#予測
trainPredict = model.predict(X_train)
testPredict = model.predict(X_test)
#print(trainPredict)
#print(testPredict)

pre = np.append(trainPredict, testPredict) #予測結果を結合
print(pre)

データのプロット

実際の平均気温変化分と予測気温変化分をそれぞれ同一グラフに表示を行いました。

#グラフ作成
plt.plot(ave_temp["年"], ave_temp["平均気温"].diff(), color = 'blue') # 平均気温(差分)をプロット
plt.plot(ave_temp["年"], pre, color = 'red') # 予測気温をプロット

plt.title('Ridge Line') # 図のタイトル
plt.xlabel('t') # x軸のラベル
plt.ylabel('Average temp')  # y軸のラベル
plt.grid()  # グリッド線を表示

plt.show()   # 図の表示

X軸が年数のみの反映でグラフが簡略化されているため、年と月で表示し直します。

#年月データとして読み込み
t = pd.read_csv("yo_tempe_data.csv",encoding = 'shift_jis',usecols=[0])
t = t.drop(range(0,120))
t = t.drop(range(402,418))

#グラフ作成
plt.plot(t["年月"], ave_temp["平均気温"].diff(), color = 'blue') # 平均気温(差分)をプロット
plt.plot(t["年月"], pre, color = 'red') # 予測気温をプロット

plt.title('Ridge Line') # 図のタイトル
plt.xlabel('t') # x軸のラベル
plt.ylabel('Average temp')  # y軸のラベル
plt.grid()  # グリッド線を表示

plt.show()   # 図の表示

縦軸が平均気温変化分、横軸が年月（1997年1月〜2020年6月）です。青色のデータが実際の平均気温変化分で、赤色が予測平均気温変化分ですが、精度はありまり高くなさそうな結果となってしまいました。

考察

結果的に、予測するまで行うことができませんでした。うまくいかなかった理由は以下と考えます。

・データ取得を限定しすぎた　　　　　　　　　　　　　　　　　　　　　　　　　　二酸化炭素濃度から予測を行おうしましたが、気象庁データが提供している二酸化炭素濃度の測定箇所が日本では3拠点しかなかったため、同拠点で二酸化炭素濃度と平均気温を測定している与那国島に絞ってしまいました。

・特徴量データが少なすぎる　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　データ取得を限定したことにより、特徴量を追加するにもメタンガス濃度データしかなく、せっかくの多変量解析というメリットを潰してしまいました。

今回の場合ではデータを限定的にしていなければ、特徴量として人口増加・ガソリン消費量・使用電力量等のデータなども追加できたので、もう少し違ったデータが出力できたのかもしれないと考えます。