[Python編] Google Earth EngineとGoogle Colaboratoryで衛星画像から火災を可視化する方法
前回は、JavaScriptで衛星データを加工して表示する方法を紹介しましたが、Pythonの画像処理や機械学習のライブラリを用いて解析したいときもあるかもしれません。
そこで、Google Earth Engine(GEE)を使ってGoogle Colabから衛星データを取得して火災を可視化してみます。
今回作成したアプリは以下に公開しました。
画面は以下のような感じです。
2022年3月23日のイルピン上空を見ています。
赤くなっているのは火災が起きていると推定した地点で、この時期に戦闘が激化していた可能性があります。
衛星データの取得と火災の可視化
GEEを利用するにはアカウントを作成する必要があります。(Googleのアカウントも要ります。)
以下のGEEのウェブサイト右上の「Sign Up」から手続きします。
(問題なければすぐに承認されて利用できるようです。)
衛星データを取得して火災を可視化するコードは以下にアップロードしました。
アカウントが作成されていれば、以下のColab notebookを実行できます。
コードの説明
ee.ImageCollectionでイメージを取得します。
データの名前(COPERNICUS/S2_SRなど)を引数に指定
filterBounds、filterDateを使って、領域や日付でデータをフィルタリング可能。
filterを使えば任意のプロパティでフィルタリング可能。(以下は、被雲率でのフィルタリング例)
ee.ImageCollection(data_name)
.filterBounds(aoi)
.filterDate(start_date, end_date)
.filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', th_cloud)))
次に、取得したデータを加工して、火災を可視化します。
GEEで取得できるSentinel-2 Level-2Aのデータは10,000倍にスケーリングされているので、注意が必要です。
火災検知のロジックには、EO Browser編で取り扱ったものと同じAFD-S2を参考にしたものを用いました。
AFD-S2は、以下の論文でHuらが提案したものです。
詳細は、以下の記事で説明していますので、割愛します。
ここではGEEのAPIを使って衛星データを加工する例を記します。
バンドの値同士の演算:add, subtract, multiply, divide, mod, powがある。
複雑な演算は、expressionに演算表現を代入して実行可能。
gt, gte, lt, lteを使って、任意の値を閾値とした二値化が可能。
以下が一例で、multiply, expression, gtを使って先述の火災検知ロジックを実装しています。
# AFD-S2をヒントにした二直線分離(rはroughの意です)
swir2_red_ratio_exp = f'(b("{SWIR2}") / b("{RED}"))'
swir2_red_ratio = img.expression(swir2_red_ratio_exp)
swir2_red_ratio_bool = swir2_red_ratio.gt(TH_SWIR2_RED_RATIO)
swir2_bool = img.select(SWIR2).gt(TH_SWIR2 * SCALING)
r_afd_s2_bool = swir2_red_ratio_bool.multiply(swir2_bool).rename('r_afd_s2_bool')
r_afd_s2_swir2 = img.select([SWIR2]).multiply(r_afd_s2_bool).rename('r_afd_s2_swir2')
r_afd_s2_nir = img.select([NIR]).multiply(r_afd_s2_bool).rename('r_afd_s2_nir')
r_afd_s2_red = img.select([RED]).multiply(r_afd_s2_bool).rename('r_afd_s2_red')
最後に、foliumで地図を作成し、その上にレイヤーを追加します。
RGBにどのバンドの値を割り振るかと輝度のスケーリングを指定します。
r_afd_s2_bool = img.select('r_afd_s2_bool').selfMask()
center = AOI.centroid(10).coordinates().reverse().getInfo()
m = folium.Map(location=center, zoom_start=ZOOM_START)
m.add_ee_layer(img,
{'bands': ['TCI_R', 'TCI_G', 'TCI_B']},
'True Color')
m.add_ee_layer(img,
{'bands': ['B12', 'B8A', 'B4'], 'min': 0, 'max': SCALING},
'SWIR', False)
m.add_ee_layer(img,
{'bands': ['r_afd_s2_swir2', 'r_afd_s2_nir', 'r_afd_s2_red'], 'min': 0, 'max': SCALING},
'rAFD-S2', False, 0.6)
m.add_ee_layer(r_afd_s2_bool,
{'palette': 'FF3616'},
'Fire')
今回はPythonの他のライブラリを使用した発展的な解析は行いませんでしたが、ee.batch.Exportを使えば、画像やメタデータなどをGoogle Driveなどに出力できるので、あとは好きなように処理できます。
一方で、GEEにも分析用の様々なAPIが用意されており、それだけでも結構高度なことができそうです。
単に土地や日付を切り替えて見て回るのが目的なら、EO Browserの方がいいと思います。
GEEについて補足
ところで、GEEのすごいところは、衛星データをローカルにダウンロードする必要が無く、Google Cloud Platform上でそのまま扱える点だと思います。Sentinel-2の衛星データは以下の公式APIからも取得可能ですが、データのサイズが1 GB近くあってダウンロードするたびに30分くらいぼんやり待っていたので、GEEは非常にありがたいです。
また、GEEのデータカタログには、Sentinel以外にも様々な衛星が登録されています。
これらを同一のAPIで取得できるのも素晴らしいです。
GEEはJavaScriptのAPIと、ブラウザ上で動くCode Editorも提供しています。
EO browserのように、結果を地図上に表示しながら開発できるので便利です。
レイヤーを切り替えたり地点や領域を登録して読み込むUIも付いていて、Pythonより楽かもしれません。
作ったアプリもそのままデプロイできるようです。
GEEのAPIリファレンスは以下です。
チュートリアルも充実しています。