[レポート] SageMaker と Autogluon と ASDI オープンデータを使用して大気汚染を予測するワークショップに参加しました / Using SageMaker, Autogluon, and ASDI open data to predict air quality #reinvent #SUS202

SageMaker や Amazon が提供する科学的オープンデータを使って大気汚染を予測するというワークショップに参加しました。
2022.12.25

この記事は公開されてから1年以上経過しています。情報が古い可能性がありますので、ご注意ください。

はじめに

こんにちは、外がマイナス 10 度でも抹茶フラペチーノが食べたくなる Classmethod Canada の Funa です。今回は、タイトルを見てなんだか面白そう!と思って参加を決めた SageMaker などを使って大気汚染を予測しようというワークショップについてレポートします。

ワークショップ概要

Using SageMaker, Autogluon, and ASDI open data to predict air quality

According to the World Health Organization (WHO), ambient air pollution leads to an estimated 4.2 million deaths per year, and around 91 percent of the world’s population lives in places where air quality levels exceed WHO limits. In this workshop, investigate the relationship between air quality and weather using Amazon SageMaker Studio. Access open data from the Amazon Sustainability Data Initiative (ASDI) with Amazon S3 and an API to explore air quality by geography. Deploy a machine learning model that uses AutoGluon (AutoML from AWS) binary classification models to predict how weather features may result in unhealthy air quality.

日本語化したものはこちらです。

世界保健機関 (WHO) によると、大気汚染により年間推定 420 万人が死亡しており、世界人口の約 91% が大気質レベルが WHO の制限を超える場所に住んでいます。 このワークショップでは、Amazon SageMaker Studio を使用して空気の質と天気の関係を調査します。 Amazon S3 と API を使用して Amazon Sustainability Data Initiative (ASDI) のオープンデータにアクセスし、地域ごとに空気の質を調査します。そして AutoGluon (AWS の AutoML) バイナリ分類モデルを使用する機械学習モデルをデプロイして、天候の特徴がどのように健康に好ましくない空気の質をもたらす可能性があるかを予測します。

ワークショップ資料

セッション内容

サステナビリティ、空気の質と天候について

Sustainability

サステナビリティ(sustainability)とは

  • 将来の世代のニーズを満たす能力を損なうことなく、現在のニーズを満たすこと(国連ブルントランド委員会、1987 年)

WHO は大気汚染を健康に対する最大の環境リスクの 1 つであると考えており、関連するリスク要因に対処することは公衆衛生を保護するための鍵となるとしています。

大気汚染物質について

公衆衛生上の懸念の最も強い空気中の汚染物質は以下であり、これらの汚染物質に短期的および長期的にさらされると健康上の問題が発生する可能性があります。

  • 粒子状物質 (PM)
    • 粗大粒子状物質 (PM10) - 粒子の直径が 0.01 mm 以下のもの
      (花粉、農業用地、採掘作業からの風で飛ばされた粉塵など)

    • 微粒子状物質 (PM2.5) - 粒子の直径が 0.0025 mm 以下のもの
      (発電施設、産業、または車両での燃料の燃焼やガス間の化学反応など)

  • 地上オゾン (O3)
  • 二酸化窒素 (NO2)
  • 二酸化硫黄 (SO2)

天気と空気の質の関連性

風、気圧、気温、日光や湿度などの天気を作り出す大気条件はオゾンと粒子状物質などの量に影響を与えるため、天気は空気の質にも影響を与えている。

Amazon Sustainability Data Initiative (ASDI) について

ASDI

ASDI とは

Amazon Sustainability Data Initiative (ASDI) では、世界中の研究者や科学者、イノベーターがサステナビリティ関連の研究を進めるのを支援するために重要な科学的データを無料で一般公開しています。

  • 現在、気象観測、海水温、気候予測データ、衛星画像など、数ペタバイトのデータをホストしている。
  • 世界中のどこにいても、ローカルストレージやコンピューティング容量に関係なく、膨大な量のデータを数分で分析できる。

今回のワークショップで使用する ASDI のデータソース

  • 天気NOAA(米国海洋大気庁) の Global Surface Summary of the Day (GSOD)
    • 1 時間ごとの詳細な観測から得られた 18 の地表気象要素に関する 1 日のデータの要約。
    • Amazon S3 上で個々のオブジェクトとしてアクセスされる。
  • 空気の質OpenAQ
    • 指定された場所 (緯度、経度) 付近の要求された大気汚染物質の 1 日平均のデータ。
    • Amazon S3 でオブジェクトとしてアクセスされるが、このワークショップでは HTTPS API 経由でアクセスを行う。

空気の質を予測するモデルを構築しよう

使用するものは以下の通りです。

  • Amazon SageMaker Studio: 機械学習用の統合開発環境 (IDE)
  • ASDI: OpenAQ や NOAA の天気のデータといったオープンデータソース
  • AutoGluon: 表形式のデータ用の AutoML フレームワーク

ワークショップ内容

家に帰ってもう一度やってみた

re:Invent でこれが初めてのワークショップ参加でかなり雰囲気に飲まれてしまったので、家でもう一度やってみました。

Amazon Sagemaker Studio ドメインとユーザーを作成する

  • ここ の記載に沿って SageMaker Studio を始める準備を行います。

Notebook を作成する

  • SageMaker Studio を立ち上げたら、File > New > Notebook を選択し、イメージとして PyTorch 1.12 Python 3.8(optimized for CPU)、カーネルとして Python3 を選択して Select をクリックします。

Notebook の作成

  • Notebook が作成されたら、Notebook 名のところにカーソルを合わせて右クリックをして air_quality_by_weather.ipynb に名前を変更します。

インストール済みの Python パッケージの確認

  • %pip list でインストールされている Python パッケージの一覧を確認します。

(結果)

Package                    Version
-------------------------- ---------------
asttokens                  2.0.8
attrs                      21.4.0
awscli                     1.25.85
(省略)

Autogluon と ipywidgets をインストールする

  • 以下のコマンドで表形式の予測を行うための Autogluon と HTML ウィジェットの ipywidgets をインストールします。
%pip install --quiet ipywidgets==7.0

%pip install --quiet autogluon

%pip show autogluon

%pip show ipywidgets

(結果) Autogluon と Ipywidgets をインストール

Python 依存関係をインポートする

  • 以下のコードを使用して Python 依存関係をインポートします。
# Import statements for packages used...
import glob
import json
import os
import shutil
import sys
from datetime import datetime
from io import StringIO
from types import SimpleNamespace

import boto3
import ipywidgets as widgets
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import requests
from autogluon.tabular import TabularPredictor
from botocore import UNSIGNED
from botocore.config import Config
from IPython.display import clear_output
from sklearn.metrics import ConfusionMatrixDisplay, confusion_matrix
from sklearn.model_selection import train_test_split

# The following is required for matplotlib plots to display in some envs...
%matplotlib inline

結果は表示されません。

クラスと変数を作成する

  • ここ の記載に沿って AQParam、AQScenario、AQbyWeatherApp クラスを作成します。

(結果)

Classes and Variables are ready.

バンクーバーのシナリオでやってみた

空気質パラメータと空気質シナリオを確認する

  • 空気質パラメータと空気質シナリオを確認します。

調査したい場所と空気質パラメータの組み合わせをシナリオと呼ぶそうです。
ここで私は独自シナリオとしてバンクーバーのシナリオなどを以下のように入れ込んで実行しました。

# Review the pre-defined AQParams and AQScenarios in this cell. You can edit these and/or use your own...
# AQParams are added with default thresholds, which can be overridden on a per-AQScenario basis.
# These AQParams are based on the OpenAQ /parameters API call where isCore=true (https://api.openaq.org/v2/parameters).
# Default thresholds where provided using data from EPA.gov (https://www.epa.gov/criteria-air-pollutants/naaqs-table).
# Confirm and adjust params or thresholds as needed for your needs... Not for scientific or health purposes.

# Instantiate main App class with explicit mlTargetLabel and mlEvalMetric provided...
AQbyWeather = AQbyWeatherApp(mlTargetLabel="isUnhealthy", mlEvalMetric="accuracy")

# Define and add new AQParams...
AQbyWeather.addAQParam(
    AQParam(1, "pm10", "µg/m³", 150.0, "Particulate Matter < 10 micrometers")
)
AQbyWeather.addAQParam(
    AQParam(2, "pm25", "µg/m³", 12.0, "Particulate Matter < 2.5 micrometers")
)
AQbyWeather.addAQParam(AQParam(7, "no2", "ppm", 100.0, "Nitrogen Dioxide"))
AQbyWeather.addAQParam(AQParam(8, "co", "ppm", 9.0, "Carbon Monoxide"))
AQbyWeather.addAQParam(AQParam(9, "so2", "ppm", 75.0, "Sulfur Dioxide"))
AQbyWeather.addAQParam(AQParam(10, "o3", "ppm", 0.070, "Ground Level Ozone"))

# Define available AQ Scenarios for certain locations with their associated NOAA GSOD StationID values...
# NOAA GSOD Station Search: https://www.ncei.noaa.gov/access/search/data-search/global-summary-of-the-day
# NOTE: For Ozone Scenarios, we're generally using 0.035 ppm to override the default threshold.
AQbyWeather.addAQScenario(
    AQScenario("vancouver", "71892099999", AQbyWeather.aqParams["pm25"], None)
)
AQbyWeather.addAQScenario(
    AQScenario("vancouver", "71892099999", AQbyWeather.aqParams["pm10"], None)
)  # Attempt at pm10 prediction.
AQbyWeather.addAQScenario(
    AQScenario("vancouver", "71892099999", AQbyWeather.aqParams["o3"], 0.035)
)
AQbyWeather.addAQScenario(
    AQScenario("toronto", "71265099999", AQbyWeather.aqParams["pm25"], None)
)
AQbyWeather.addAQScenario(
    AQScenario("toronto", "71265099999", AQbyWeather.aqParams["o3"], 0.035)
)
AQbyWeather.addAQScenario(
    AQScenario("edmonton", "71879099999", AQbyWeather.aqParams["pm25"], None)
)
AQbyWeather.addAQScenario(
    AQScenario("edmonton", "71879099999", AQbyWeather.aqParams["o3"], 0.035)
)
AQbyWeather.addAQScenario(
    AQScenario("winnipeg", "71852099999", AQbyWeather.aqParams["pm25"], None)
)
AQbyWeather.addAQScenario(
    AQScenario("winnipeg", "71852099999", AQbyWeather.aqParams["o3"], 0.035)
)
AQbyWeather.addAQScenario(
    AQScenario("montreal", "71612099999", AQbyWeather.aqParams["pm25"], None)
)
AQbyWeather.addAQScenario(
    AQScenario("montreal", "71612099999", AQbyWeather.aqParams["o3"], 0.035)
)
AQbyWeather.addAQScenario(
    AQScenario("calgary", "71877099999", AQbyWeather.aqParams["pm25"], None)
)
AQbyWeather.addAQScenario(
    AQScenario("calgary", "71877099999", AQbyWeather.aqParams["o3"], 0.035)
)
AQbyWeather.addAQScenario(
    AQScenario("saskatoon", "71496099999", AQbyWeather.aqParams["pm25"], None)
)
AQbyWeather.addAQScenario(
    AQScenario("halifax", "71395099999", AQbyWeather.aqParams["pm25"], None)
)

print(f"AQbyWeather.aqParams: {str(len(AQbyWeather.aqParams))}")
print(
    f"AQbyWeather.aqScenarios: {str(len(AQbyWeather.aqScenarios))} (Default Selected: {AQbyWeather.selectedScenario.name})"
)

(結果)

AQbyWeather.aqParams: 6
AQbyWeather.aqScenarios: 15 (Default Selected: vancouver_pm25)

モデル化するシナリオを選択する

  • ここで以下を実行して vancouver_pm25 シナリオを選択しました。
# Select a Scenario via DROP DOWN LIST to use throughout the Notebook. This will drive the ML process...
# A default "value" is set to avoid issues. Change this default to run the Notebook from start-to-finish for that Scenario.
print("*** CHOOSE YOUR OWN ADVENTURE HERE ***")
print("Please select a Scenario via the following drop-down-list...")
print("(NOTE: If you change Scenario, you must re-run remaining cells to see changes.)")
ddl = widgets.Dropdown(
    options=AQbyWeather.aqScenarios.keys(),
    value=AQbyWeather.aqScenarios["vancouver_pm25"].name,
)  # <-- DEFAULT / FULL-RUN VALUE
ddl

(結果)

*** CHOOSE YOUR OWN ADVENTURE HERE ***

Please select a Scenario via the following drop-down-list...

(NOTE: If you change Scenario, you must re-run remaining cells to see changes.)

シナリオを実行する

  • 以下のコードを使用して先ほど選択したシナリオを実行します。

これで ML モデルが使用するターゲットデータが入力されるそうです。

if ddl.value:
    AQbyWeather.selectedScenario = AQbyWeather.aqScenarios[ddl.value]
    print(AQbyWeather.selectedScenario.getSummary())
    print(AQbyWeather.selectedScenario.toJSON())
else:
    print("Please select a Scenario via the above drop-down-list.")

(結果)

Scenario: vancouver_pm25 => Particulate Matter < 2.5 micrometers (pm25) with UnhealthyThreshold > 12.0 µg/m³
{
  "aqParamTarget": {
    "desc": "Particulate Matter < 2.5 micrometers",
    "id": 2,
    "name": "pm25",
    "unhealthyThresholdDefault": 12.0,
    "unit": "\u00b5g/m\u00b3"
  },
  "aqRadiusMeters": 16100,
  "aqRadiusMiles": 10,
  "location": "vancouver",
  "modelFolder": "AutogluonModels",
  "name": "vancouver_pm25",
  "noaaStationID": "71892099999",
  "noaaStationLat": 0.0,
  "noaaStationLng": 0.0,
  "openAqLocIDs": [],
  "unhealthyThreshold": 12.0,
  "yearEnd": 2022,
  "yearStart": 2016
}

NOAA 気象データを取得する

  • 以下のコードを使用して NOAA 気象データを取得します。
# GET NOAA GSOD WEATHER DATA...
print(AQbyWeather.selectedScenario.getSummary())
noaagsod_df = AQbyWeather.getNoaaDataFrame()
noaagsod_df = noaagsod_df[~noaagsod_df.eq(9999.9).any(1)]

if len(noaagsod_df) >= 1:
    # Update NOAA Station Lat/Lng...
    AQbyWeather.selectedScenario.updateNoaaStationLatLng(noaagsod_df.iloc[0])

    # Save DataFrame to CSV...
    noaagsod_df.to_csv(AQbyWeather.getFilenameNOAA(), index=False)

    # Output DataFrame properties...
    print("noaagsod_df.shape =", noaagsod_df.shape)
    display(noaagsod_df)
else:
    print("No data is available for this location.")

すると、2016 年 1 月から 2022 年 12 月までのバンクーバーの気象データが返却されました。
(ちょうど 2016 年 1 月末からバンクーバーに住んでいるので、なんだか感慨深い気持ちになりました。)

2016 年 1 月から 2022 年 12 月までのバンクーバーの気象データ

※ 気象データのラベルの意味はこちらで確認できます。

OpenAQ 空気質の日毎の平均データを取得する

  • 以下のコードを使用して OpenAQ 空気質の日毎の平均データを取得します。
# GET OPENAQ AIR QUALITY DAILY AVERAGES DATA...
print(AQbyWeather.selectedScenario.getSummary())
aq_df = (
    AQbyWeather.getOpenAqDataFrame()
)  # Gets nearby Location IDs THEN gets associated daily averages.

if len(aq_df) > 0:
    # Output DataFrame properties...
    print("aq_df.shape =", aq_df.shape)
    display(aq_df)
    aq_df.to_csv(AQbyWeather.getFilenameOpenAQ(), index=False)
else:
    print("No data is available for this location.")

そうすると、2016 年 1 月から 2022 年 12 月までのバンクーバーの空気質の日毎の平均データが返却されました。
(2018 年の 8 月がめちゃくちゃ unhealthy な結果となっています...あの頃何してただろう...)

2016 年 1 月から 2022 年 12 月までのバンクーバーの空気質の日毎の平均データ

NOAA と OpenAQ のデータを 1 つのデータフレームにマージする

  • 以下のコードを使用して NOAA と OpenAQ のデータを 1 つのデータフレームにマージします。

isUnhealthy が 1 の場合は大気の質が不健康だったことを表し、0 は大気の質が健全だったことを表します。

# Merge the NOAA GSOD weather data with our OpenAQ data by DATE...
# Perform another column drop to remove columns we don't want as features/inputs.
# This column removal will NOT be necessary once we can use Autogluon ignore_columns param (TBD).
print(AQbyWeather.selectedScenario.getSummary())
merged_df = AQbyWeather.getMergedDataFrame(noaagsod_df, aq_df)

if len(merged_df > 0):
    # Output DataFrame properties...
    print("merged_df.shape =", merged_df.shape)
    display(merged_df)
    merged_df.groupby([AQbyWeather.mlTargetLabel]).size().plot(kind="bar")
    merged_df.to_csv(AQbyWeather.getFilenameOther("dataMERGED"), index=False)

その結果、ほぼ healthy であることがわかりました(偏りすぎて確認結果としては良くないかもしれません...)。

ほぼ healthy な結果に

マージされたデータフレームで相関関係を視覚化する

  • 以下のコードを実行してマージされたデータフレームで相関関係を視覚化することで、空気の質を予測するためのモデルでどの変数が役立つかを確認します。
# Visualize correlations in our merged dataframe...
print(AQbyWeather.selectedScenario.getSummary())
correlations = merged_df.corr()
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111)
cax = ax.matshow(correlations, vmin=-1, vmax=1,cmap = "coolwarm")
fig.colorbar(cax)
ticks = np.arange(0, len(merged_df.columns), 1)
ax.set_xticks(ticks)
ax.set_yticks(ticks)
ax.set_xticklabels(merged_df.columns)
ax.set_yticklabels(merged_df.columns)
plt.show()

相関プロットは以下の通りです。この図を見ると、AirQuality (isUnhealthy) は WindSpeed とわずかに負の相関があり、Month とわずかに相関しています。これは、通常、風速が高いほど空気の質が高く、空気の質が季節に関連していることを意味しているそうです...!

相関プロット

最後までやったけど...

本当はこのワークショップにはまだ続きがあって、家に帰ってから最後までやったのですが、内容がとても難しくて今の自分にはかなり厳しかったです。機械学習についての知識をつけてから出直してきます。

感想

バンクーバーは以下の傾向があります。

  • 雨季がやたら長く、冬は他のカナダの都市と比べて寒くなりすぎない。
  • あっという間に過ぎ去るすごしやすい夏、自然が多いので山火事が発生する。
  • 地価が激しく高いためか、(おそらく)産業地区がない、環境対策に力を入れている。

夏に山火事が起きた時くらいしか、外にいて空気の質が悪いなと思うことはないので空気の質は普通〜割と良いのではないかと思っていましたが、それは当たっていたようです。

自分の生活に密接に関わることが機械学習で予測できるというのは非常に面白くて、普段は全く機械学習とは接点がないのですが、今後もちょくちょく手を出せたら良いなと思いました。