【Python】横浜市の花粉飛散量予測モデル

0.目次

1．はじめに
2．開発環境
3．使用したデータ
4．予測モデル作成内容
　4－1．データの取得
　4－2．特徴量の追加
　4－3．訓練用データと検証用データ
　4－4．アルゴリズム
　4－5．評価
5．コード
　5－1．インポート、データのダウンロード
　5－2．データの前処理
　5－3．特徴量の追加
　5－4．特徴量の分析
　5－5．訓練データと検証用データに分割
　5－6．LightGBMで予測、評価
6．予測結果
7．おわりに
8．参考資料

1．はじめに

Aidemyのデータ分析講座の最終課題として毎春自分が悩まされている花粉の飛散量を予測するモデルを作成してみました。

2．開発環境

Windows 11 Home
Google Colaboratory
Python 3.9.16

3．使用したデータ

環境庁花粉観測システムの花粉観測データ

• 期間：2012年～2020年（2月～6月）
• 拠点：横浜市
• データ内容：1時間ごとの花粉飛散量、降雨量、気温、風向、風速

気象庁ホームページの気象データ

• 期間：2012年～2020年（7月～8月）
• 拠点：横浜市
• データ内容：1時間ごとの降雨量、気温

4．予測モデル作成内容

4－1．データの取得

環境庁花粉観測システムのデータには2月から6月のデータのみが含まれています。(6月は花粉データは全てゼロになっています。)一般的に花粉の飛散量は前年の夏の気温や降水量と関係していると言われているので、気象庁のホームページから同じ期間の7月、8月の気温と降水量のデータをダウンロードして追加しました。

4－2．特徴量の追加

以下を特徴量として定義し、Seabornのヒートマップでデータ同士の相関関係を確認したり、実際に何度もモデル作成を行って予測結果を確認し、最終的に【使用】と書いた項目のみを使用しました。
　・前年2月から5月の積算花粉飛散量【使用】
　・前年6月から8月の積算気温【使用】
　・前年6月から8月の積算降水量【使用】
　・直近1週間の平均花粉飛散量【使用】
　・直近1週間の平均気温
　・直近1週間の平均降水量【使用】
　・直近1日の平均花粉飛散量【使用】
　・直近1日の平均気温
　・直近1日の平均降水量【使用】
　・直近12時間の平均花粉飛散量【使用】
　・直近12時間の平均気温
　・直近12時間の平均降水量【使用】
　・直近6時間の平均花粉飛散量【使用】
　・直近6時間の平均気温【使用】
　・直近6時間の平均降水量【使用】

ヒートマップを作成したところ、以下のように今回予測したい対象である花粉(pollen)とその他の項目の相関関係が表示されました。(項目同士の関係が深いほど薄い色で表示されます。)

4－3．訓練用データと検証用データ

• 2013年から2019年のデータを訓練用データとして使用
  （2012年のデータには前年の夏の気温や降水量データが含まれていないため）
• 2020年3月1日から3月15日のデータを検証用データとして使用

4－4．アルゴリズム

今回はLightGBM(regression_L1)を使用してモデリングを行いました。
LightGBMは、勾配ブースティングと呼ばれる手法を用いたアルゴリズムの一種で、RandomForestと同様、決定木がもとになっています。LightGBMでは、前に作った決定木の誤差が小さくなるように次の決定木を作成します。RandomForestが並列ならば、LightGBMは直列に決定木を作成するイメージです。

＜LightGBMイメージ図＞※Aidemy教材より引用

花粉の飛散量のデータは、複雑な形をしているので、複雑さに対応できそうなモデルであること、またデータ量も多い（45,865行）ので、短時間でモデリングできることからLightGBMを採用しました。（実際には他にも複数のアルゴリズムを試しましたが、データの精度が非常に低かったり、時間がかかりすぎてなかなか結果を得ることができなかっため、LightGBMのみを使用しました。）

4－5．評価

今回はKFoldでクロスバリデーションを行いました。KFoldは、用意したデータセットをk個に分割し、 そのうちの1つを評価用データ、残りのk-1個を訓練データとして使用します。分割した分だけ、学習と評価を繰り返し、得られるk個のモデルと性能評価から平均性能を算出します。この手法は、「複数回データを変えて実行するので安定した評価を行うことができる」、「全てのデータを学習に用いるのでデータを効率よく利用し汎化性能の高いモデルが作成できる」などの利点があります。
※scikitLearnドキュメントから引用

また、評価指標としてMean Absolute Error（絶対平均誤差）を使用しています。
Mean Absolute Errorは真の値と予測値の差の絶対値の平均です。花粉量のデータは時々大量に飛散する外れ値があるので、外れ値に影響されにくい評価指標を選びました。

5．コード

5－1．インポート、データのダウンロード

ライブラリをインポートし、あらかじめGoogleドライブに保存しておいたデータをダウンロードした。

インポート、データのダウンロード

#import
from google.colab import drive
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%pip install seaborn
import seaborn as sns
from sklearn.model_selection import KFold
import lightgbm as lgb
from sklearn.metrics import mean_absolute_error

#データの読み込み
from google.colab import drive
drive.mount('/content/drive/')
data = pd.read_csv("/content/drive/MyDrive/Data/pollen_yokohama_2012-2020_v4.csv")

5－2．データの前処理

データの欠測値を0に置き換えたり、文字列を数値に置き換える処理を行った。

データの前処理

#データの前処理
#花粉飛散量がマイナスのもの（黄砂や降雪などによる欠測）を0に置換
data.loc[lambda df: df["pollen"] < 0, "pollen"] = 0

#花粉飛散量をグラフで表示
#data['pollen'].plot(grid=True)

#"欠測"(気温、風量等)をNaNに置換
data = data.replace('欠測', 0)

# 気温、風向、風速の値を文字列から数値へ変更
data = data.astype({'temparature' : float,'winddirection' : int,'windspeed' : float})

5－3．特徴量の追加

特徴量の追加を行った。

#昨年の2月から5月の積算花粉飛散量

sum_pollen_2012= data[(data["year"]==2012) & (data["month"]>=2) & (data["month"]<=5)].sum(axis=0)["pollen"]
sum_pollen_2013= data[(data["year"]==2013) & (data["month"]>=2) & (data["month"]<=5)].sum(axis=0)["pollen"]
sum_pollen_2014= data[(data["year"]==2014) & (data["month"]>=2) & (data["month"]<=5)].sum(axis=0)["pollen"]
sum_pollen_2015= data[(data["year"]==2015) & (data["month"]>=2) & (data["month"]<=5)].sum(axis=0)["pollen"]
sum_pollen_2016= data[(data["year"]==2016) & (data["month"]>=2) & (data["month"]<=5)].sum(axis=0)["pollen"]
sum_pollen_2017= data[(data["year"]==2017) & (data["month"]>=2) & (data["month"]<=5)].sum(axis=0)["pollen"]
sum_pollen_2018= data[(data["year"]==2018) & (data["month"]>=2) & (data["month"]<=5)].sum(axis=0)["pollen"]
sum_pollen_2019= data[(data["year"]==2019) & (data["month"]>=2) & (data["month"]<=5)].sum(axis=0)["pollen"]
data.loc[data['year'] == 2013, 'sum_pollen_ly'] = sum_pollen_2012
data.loc[data['year'] == 2014, 'sum_pollen_ly'] = sum_pollen_2013
data.loc[data['year'] == 2015, 'sum_pollen_ly'] = sum_pollen_2014
data.loc[data['year'] == 2016, 'sum_pollen_ly'] = sum_pollen_2015
data.loc[data['year'] == 2017, 'sum_pollen_ly'] = sum_pollen_2016
data.loc[data['year'] == 2018, 'sum_pollen_ly'] = sum_pollen_2017
data.loc[data['year'] == 2019, 'sum_pollen_ly'] = sum_pollen_2018
data.loc[data['year'] == 2020, 'sum_pollen_ly'] = sum_pollen_2018

#昨年の6月から8月積算気温

sum_temparature_2012 =data[(data["year"]==2012) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["temparature"]
sum_temparature_2013 =data[(data["year"]==2013) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["temparature"]
sum_temparature_2014 =data[(data["year"]==2014) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["temparature"]
sum_temparature_2015 =data[(data["year"]==2015) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["temparature"]
sum_temparature_2016 =data[(data["year"]==2016) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["temparature"]
sum_temparature_2017 =data[(data["year"]==2017) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["temparature"]
sum_temparature_2018 =data[(data["year"]==2018) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["temparature"]
sum_temparature_2019 =data[(data["year"]==2019) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["temparature"]
data.loc[data['year'] == 2013, 'sum_temparature_ly'] = sum_temparature_2012
data.loc[data['year'] == 2014, 'sum_temparature_ly'] = sum_temparature_2013
data.loc[data['year'] == 2016, 'sum_temparature_ly'] = sum_temparature_2015
data.loc[data['year'] == 2017, 'sum_temparature_ly'] = sum_temparature_2016
data.loc[data['year'] == 2018, 'sum_temparature_ly'] = sum_temparature_2017
data.loc[data['year'] == 2019, 'sum_temparature_ly'] = sum_temparature_2018
data.loc[data['year'] == 2020, 'sum_temparature_ly'] = sum_temparature_2019

#昨年の6月から8月の積算降水量

sum_precipitaion_2012 =data[(data["year"]==2012) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["precipitation"]
sum_precipitaion_2013 =data[(data["year"]==2013) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["precipitation"]
sum_precipitaion_2014 =data[(data["year"]==2014) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["precipitation"]
sum_precipitaion_2015 =data[(data["year"]==2015) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["precipitation"]
sum_precipitaion_2016 =data[(data["year"]==2016) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["precipitation"]
sum_precipitaion_2017 =data[(data["year"]==2017) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["precipitation"]
sum_precipitaion_2018 =data[(data["year"]==2018) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["precipitation"]
sum_precipitaion_2019 =data[(data["year"]==2019) & (data["month"]>=6) & (data["month"]<=8)].sum(axis=0)["precipitation"]
data.loc[data['year'] == 2013, 'sum_precipitation_ly'] = sum_precipitaion_2012
data.loc[data['year'] == 2014, 'sum_precipitation_ly'] = sum_precipitaion_2013
data.loc[data['year'] == 2015, 'sum_precipitation_ly'] = sum_precipitaion_2014
data.loc[data['year'] == 2016, 'sum_precipitation_ly'] = sum_precipitaion_2015
data.loc[data['year'] == 2017, 'sum_precipitation_ly'] = sum_precipitaion_2016
data.loc[data['year'] == 2018, 'sum_precipitation_ly'] = sum_precipitaion_2017
data.loc[data['year'] == 2019, 'sum_precipitation_ly'] = sum_precipitaion_2018
data.loc[data['year'] == 2020, 'sum_precipitation_ly'] = sum_precipitaion_2019

#直近一週間の平均花粉量、平均気温、平均降水量

#直近1週間の平均花粉飛散量
data["ave_pollen_1w"] = data["pollen"].rolling(168).mean()

#直近1週間の平均気温
data["ave_temparature_1w"]=data["temparature"].rolling(168).mean()

#直近1週間の平均降水量
data["ave_pricipitation_1w"] =data["precipitation"].rolling(168).mean()

#直近24時間の平均花粉量、平均気温、平均降水量

#直近24時間の平均花粉飛散量
data["ave_pollen_1d"] = data["pollen"].rolling(24).mean()

#直近24時間の平均気温
data["ave_temparature_1d"]=data["temparature"].rolling(24).mean()

#直近24時間の平均降水量
data["ave_pricipitation_1d"] =data["precipitation"].rolling(24).mean()

#直近12時間の平均花粉量、平均気温、平均降水量

#直近12時間の平均花粉飛散量
data["ave_pollen_12h"] = data["pollen"].rolling(12).mean()

#直近12時間の平均気温
data["ave_temparature_12h"]=data["temparature"].rolling(12).mean()

#直近12時間の平均降水量
data["ave_pricipitation_12h"] =data["precipitation"].rolling(12).mean()

#直近6時間の平均花粉量、平均気温、平均降水量

#直近6時間の平均花粉飛散量
data["ave_pollen_6h"] = data["pollen"].rolling(6).mean()

#直近6時間の平均気温
data["ave_temparature_6h"]=data["temparature"].rolling(6).mean()

#直近6時間の平均降水量
data["ave_pricipitation_6h"] =data["precipitation"].rolling(6).mean()

5－4．特徴量の分析

項目同士の相関関係を確認するためにヒートマップを表示。

ヒートマップを表示

plt.figure(figsize=(13, 11))
sns.heatmap(data.corr())
plt.tight_layout();

5－5．訓練用データと検証用データに分割

訓練用データと検証用データに分割し、訓練に使用しないデータを削除した。

データの分割

#6月から8月のデータには花粉の飛散量がないので訓練データから削除
data=data[data['month'] != 6]
data=data[data['month'] != 7]
data=data[data['month'] != 8]

#2012年のデータには昨年度の積算花粉飛散量等がないので訓練データから削除
data=data[data['year'] != 2012]

#データを訓練用データと検証用データに分割
train = data[data['testflag']==0]
test = data[data['testflag']==1]
train_target = train['pollen']
test_target = test['pollen']

#今回訓練に用いないを項目を削除'
drop_col=['month','pollen','testflag','ave_temparature_1w','ave_temparature_1d','ave_temparature_12h']

train = train.drop(drop_col, axis=1)
test = test.drop(drop_col, axis=1)

5－6．LightGBMで予測、評価

LightGBM

#LightGBMで予測

#K-Foldを定義
cv = KFold(n_splits=3, random_state=0, shuffle=True)

train_acc_list = []
val_acc_list = []
models = []

# LightGBMのパラメータを定義　
lgb_params = {
    "objective":"regression_l1",
    "metric":"mean_absolute_error",
    "force_row_wise" : True,
    "seed" : 0,
    'learning_rate': 0.1,
    'min_data_in_leaf': 5
    }

#モデルを学習させ、交差検証法で評価
for i ,(trn_index, val_index) in enumerate(cv.split(train, train_target)):
    
    print(f'Fold : {i}')
    X_train ,X_val = train.iloc[trn_index], train.iloc[val_index]
    y_train ,y_val = train_target.iloc[trn_index], train_target.iloc[val_index]
    
    # LigthGBM用のデータセットを定義
    lgb_train = lgb.Dataset(X_train, y_train)
    lgb_valid = lgb.Dataset(X_val, y_val)
    
    model = lgb.train(
        params = lgb_params, 
        train_set = lgb_train,
        valid_sets = [lgb_train, lgb_valid], 
        verbose_eval = 100 ,
        early_stopping_rounds=10
       )
    
    y_pred = model.predict(X_train)
    train_acc = mean_absolute_error(y_train,y_pred)
    train_acc_list.append(train_acc)
    
    y_pred_val = model.predict(X_val)
    val_acc = mean_absolute_error(y_val,y_pred_val)
    val_acc_list.append(val_acc)
    
    models.append(model)

print('-'*10 + 'Result' +'-'*10)
print(f'Train_acc : {train_acc_list} , Ave : {np.mean(train_acc_list)}')
print(f'Valid_acc : {val_acc_list} , Ave : {np.mean(val_acc_list)}')


#2020年の飛散量の予測
test_pred = model.predict(test)
test_acc =mean_absolute_error(test_target, test_pred)

print(test_acc)

#グラフ作成準備
pred =pd.DataFrame(test_pred,columns=['pred'])#予測値をDataFrameに格納

test["test_target"]= test_target
test = test.reset_index()
test_result=pd.concat([test,pred], axis=1)

#予測値および正解の値（実測値）をグラフに出力
test_result['test_target'].plot(grid=True,color ="b")
test_result['pred'].plot(grid=True,color ="r")
plt.show()

6．予測結果

結果は以下の通りでした。

• 訓練用データの絶対平均誤差（の平均）：14.26
• バリデーションデータの絶対平均誤差（の平均）：15.60
• 検証用データの絶対平均誤差：9.89

絶対平均誤差の数値だけではわかりにくいので、検証用のデータをグラフにしてみました。

• 実測値（正解）：青
• 予測値：赤

花粉量が突出している部分がどれだけ突出しているかについては予測できていませんが、大まかな花粉量の推移や、花粉量が大量飛散するタイミングについてはほぼ予測することができました。

花粉の飛散量は現在は環境庁では公開していませんが、ウェザーニュース社のサイトでほぼリアルタイムのデータを取得することができますし、気象庁のサイトで気温や降水量などのデータもほぼリアルタイムで取得することができるので、これらを活用すれば花粉症を患う人にとって十分有益な結果が得られると思います。

7．おわりに

私はPythonの初心者で、本格的に機械学習のコードを書くのは今回が初めてだったので、基本的な部分（グラフを表示させるためにデータを統合する等）に時間がかかりましたが、モデリング自体はほぼ自動的にプログラムが行ってくれるので、（私の習熟度から考えると）驚くほど簡単に精度の高い結果を得られたと感じています。

今回のように一時間単位での花粉の飛散量の予測を行うためには、前年のデータ（気温や降水量等）よりも、直近の花粉飛散量や降水量の方が影響が大きいことも何度もモデリングをする中で良く分かりました。（前年の飛散量や夏の気温、降水量のみを特徴量として追加した際の結果は以下でした。花粉量の増減は何となく予測されていますが、花粉が大量飛散するタイミング等については予測できていません。）

Aidemyの講師の先生にヒートマップについて教えていただいて、どの特徴量が有効か簡単に見分けられるようになったことは精度の向上に大きく役立ちました。（直近の花粉量が重要であることもヒートマップを見れば分かります。）

今回は突出して花粉の飛散量が多くなる部分が、どれくらい多くなるのかについて予測できなかったので、今後はこのようなデータをどのようにすれば予測できるか、学んでいきたいと思います。

8．参考資料

森林総合研究所　平成21年版　研究成果選集「スギ花粉はどこから飛んでくるのか？」
ウェザーニューズに聞く、花粉の飛散量予測の最前線とビジネスチャンス
花粉飛散量予測コンペティション
【Python】seabornのグラフを活用したデータ分析の手法メモ