コロナ自粛で空気がキレイになった地域を Google Earth Engine × Mapbox で可視化する
これは何の記事？（TL;DR）
コロナの自粛期間で大気汚染物質（NO2）が減った地域を可視化しました。
実装までのプロセスは、
Google Earth Engineで衛星データを取得・分析 -> Mapboxで可視化
といった流れです。
この記事では、この可視化までの道のりを紹介します。
はじめに
少し前、新型コロナウイルスの影響で大気汚染が改善していることが、様々なメディアで報道されました。
例えば、下記のBBCやEarth Observatoryの記事では、都市のロックダウンが大気汚染を改善させたとの内容が示されています。
日本でも、ロックダウンはされていないものの、緊急事態宣言による全国的な移動自粛期間がありました。
つまり、日本でも大気汚染が改善された可能性は十分にあります。
また、先日NASAらは、コロナの影響を監視するためのダッシュボードを公開しました。
このダッシュボードでは主に、次の環境を世界の一部地域で監視しています。
- 水質（クロロフィル）
- 大気（NO2、CO2）
- 夜間光量
- 空港や駐車場の混み具合
ただ、日本で情報が確認できる地域は「東京エリア」のみです。（2020/07/11 現在）
そこで、日本全域で大気がコロナ自粛でどう変わったかを、Google Earth EngineとMapboxを利用して可視化してみました。
今回は、その方法を紹介します。
何を可視化するか
日本の緊急事態宣言が大気汚染にどう影響したか確かめるべく、大気汚染物質の１つとされるNO2の前年比を可視化します。
比較期間は、緊急事態宣言の対象地域が全国となった2020年4月16日*〜１ヶ月間と、前年の同期間です。
（※比較期間は、前年だけでなく複数年を対象とするのが望ましいが、今回はデータセットの都合で前年のみとした）
対象読者
この記事は、次のような読者を想定しています。
- 衛星データの解析を手軽にやりたい人
- Google Earth Engineを使い始めた／これから使う人
- Mapboxに載せるデータの幅を広げたい人
Google Earth Engine (GEE) とは
Googleがクラウド管理する大規模な衛星データを分析・出力できるプラットフォームです。
衛星データは検索・取得・整形が大変で、分析にもかなりのマシンパワーが必要ですが、これら一連の作業をカジュアルにできるようにしたのが、このGEEです。
GEE の利用条件
GEEは、研究／教育／非営利目的での使用に関して無料で利用できます。商用利用に関しては、有料の商用ライセンスが提供されています。
GEE の利用方法
GEEの利用方法はいくつかありますが、
で利用されるのが一般的です。
Earth Engine Code Editorは、地図ビュアーと連動したUIで、分析結果を地図上に載せるインタラクティブな分析をする開発者に愛されています。
またEarth Engine Python APIは、Jupyter notebookやGoogle Colabから呼び出して、別のドメイン（画像処理やNLPなど）の計算結果と組み合わせたい研究者にニーズがあります。
ちなみにこの記事では、前者のCode Editorを使います。
手順
1. GEEで衛星データを取得・分析・出力する
まずはGEEにて衛星データを処理します。
今回はEarth Engine Code Editorを使います。
（初めて利用する際はこちらからサインアップが必要です）
1-1. 衛星データを取得
手始めに、Code EditorでGoogleのクラウド上にある衛星データにアクセスして、表示するところまでをやってみましょう。
次のコードを実行してみてください。
実行すると、Sentinel-5から観測するNO2のデータセットより、２つの期間（2019年／2020年）のデータがそれぞれ表示されます。
1-2. 取得した衛星データを分析
1-1.で取得した２つの期間の衛星データの差分を取ってみましょう。
上記コードを実行すると、2020年のNO2データが2019年と比べてどれだけ増減しているかが観測できます。
東京／大阪／愛知など、人の集まりやすい場所でNO2が減っていることが見てとれます。
1-3. 分析した衛星データを出力
1-2.の分析結果をMapboxにアップロードできる形式に変換・出力します。
今回は、次の２パターンのデータ形式に変換・出力します。
- MBTiles（→ 1-3A）
- GeoJSON（→ 1-3B）
どちらか好みの方を選んでください。
1-3A. MBTiles形式のデータを作る
MBTilesとは、正方形のPNG画像（タイル）をタイル状に配置してまとめたラスタ形式のデータです。
そして、Mapboxにアップロードできる形式の中で、容量制限が最も大きな形式です。（2020/07/11 現在）
この手順では、GEEから出力したGeoTIFFを、次の順番でMBTilesに変換します。
GeoTIFF -> Shapefile -> GeoJSON -> MBTiles
まずは、次のコードで1-2.の分析結果をGeoTIFF形式で出力しましょう。
次にこのGeoTIFFを、GDALの gdal_polygonize でShapefileに変換します。
（GDALのインストール方法はこちらの Github wiki が参考になります）
# Convert a GeoTIFF to a Shapefile $ gdal_polygonize.py input_tif_from_earth_engine.tif -mask input_tif_from_earth_engine.tif -f "ESRI Shapefile" output_shape_file output_shape_file DN
そして変換したShapefileを、GDALの ogr2ogr でGeoJSONに変換します。
# Convert a Shapefile to GeoJSONs $ ogr2ogr -f GeoJSON output_geojson_file.json -t_srs http://spatialreference.org/ref/epsg/4326/ input_shape_file.shp --config OGR_ENABLE_PARTIAL_REPROJECTION TRUE -skipfailures
これで、いくつかのスケールに対応したGeoJSON群が得られるはずです。
続いてMapboxのライブラリであるTippecanoeを使って、これらの各スケールのGeoJSONをズームレベルに割り当て、それぞれを１枚のタイルセットとしてMBTtilesに統合します。
# Create a MBTiles. (for zooms 0 through 4, on the 5,000m scale data) $ tippecanoe -z4 -o output_5k.mbtiles -pd -l output_5k input_5k.json input_5k_2.json ... # Create a MBTiles. (for zooms 5 through 6, on the 1,000m scale data) $ tippecanoe -Z5 -z6 -o output_1k.mbtiles -pd -l output_1k input_1k.json input_1k_2.json ... # Create a MBTiles. (for zooms 7 through 10, on the 250m scale data) $ tippecanoe -Z7 -z10 -o ouput_250m.mbtiles -pd -l output_250m input_250m.json input_250m_2.json input_250m_3.json ...
今回は３つのスケールでMBTilesを作成しました。
ここで、３つのMBTilesを、さらに１つのMBTilesにまとめます。
# Merge three MBTiles into one MBTile $ tile-join -o output_combined.mbtiles output_1k.mbtiles output_5k.mbtiles output_250M.mbtiles
これで完成です。
（ここまでの全コードはこちら）
1-3B. GeoJSON形式のデータを作る
GeoJSONとは、地理空間情報をJSON形式で記述するためのデータ形式で、要は型の決まったJSONです。
MBTilesはラスターデータでしたが、GeoJSONはベクトルデータです。
この手順では、GEEのCode Editorから直接GeoJSONを出力します。
ざっくり解説すると上記コードは、欠損値をマスクで取り除き、ラスター形式をベクトル形式に変換するための関数で変換したデータを、GeoJSON形式で出力するという構成です。
ただし全世界の全ピクセルを計算するには、たとえGEEでも計算パワーが足りないため、今回はデータを日本の形にクリッピングして計算量を抑えました。
こちらを実行してGeoJSONを生成できれば、この手順は完了です。
（ここまでの全コードはこちら）
2. GEEで出力した衛星データをMapboxで可視化する
2-1. GEEのデータをMapboxのカスタム地図に落とし込む
1.で作成したMBTilesまたはGeoJSONをMapboxにアップロードし、Mapbox Studioで地図上に落とし込みます。
この手順は、私が以前書いた記事に載っていますので、そちらをお読みください。
（2. Mapbox Studio でカスタム地図（Style）を作成に記載）
（デモで表示しているカスタム地図はこちら）
2-2. カスタム地図をWebで表示
2-1.にてMapbox Studioで作ったカスタム地図を、Mapbox GL JSで読み込みます。
例えば、HTMLの<head>で次のMapbox GL JSライブラリを読み込み、
<body>に次のコードを含めればOKです。
（ちなみにデモのコードはこちら）
おわりに
この記事では、Google Earth EngineとMapboxを利用して、衛星データから大気汚染物質（NO2）が減った地域を可視化する方法を紹介しました。
Google Earth EngineもMapboxも素晴らしいツールです。
Google Earth Engineには、全世界規模のデータを爆速で操作できる楽しさがあります。
今回取り上げたNO2の他にも、最新の人口統計データや高解像な衛星画像など、様々なデータセットが揃っています。
また、分類や回帰式、Tensorflowモデルの利用などにも対応しているので、汎用的な衛星分析プラットフォームとして皆さんも是非活用してみてください。
そしてMapboxは、直感操作で思い通り以上の地図が作れるツールです。
今回はあまり焦点を当てませんでしたが、実用性・デザイン性ともに最高の地図ツールです。
Twitterのハッシュタグ #BuiltWithMapbox で世界のMapboxユーザの作品が検索できるので、気になる方は是非覗いてみてください。
おまけ：可視化結果の考察
今回は単なる前年比を可視化しただけなので、強い根拠とは言えないですが、一応考察します。（デモページより）
可視化の対象としたNO2は、化石燃料を燃やした際に排出される副産物です。
自粛期間中、工場の活動や車・飛行機などの交通量が減少したと考えると、都市部でNO2が減少している結果には納得ができます。
ただ、全国的にNO2は前年より減少していますが、一部で微増している地域も見られます。
例えば、越谷市や各務原市。これらは都市部に近いベッドタウンです。
これは想像ですが、NO2が増えた原因の一つは、普段都市に出稼ぎするビジネス層がテレワークで、地元で活動する機会が増えたためではないかと考えられます。
とは言え、「他のベッドタウンは？」とも思うので、もし詳細が気になる方は次のような項目を追加調査してみてください。
- NO2以外のエアロゾル
- 車両交通量
またNO2は、多くの科学者が気候変動に寄与しているとの見方を示しています。
すぐに気候変動に直結しないとは思いますが、長い目で見て...
- 地表温度
- 海面温度
などとの相関を調べてみると、面白い結果が出るかもしれません。
