🛰️

速攻単純 Sentinel-2 光学衛星の雲除去をする方法

2023/12/13に公開

導入

光学衛星画像を取得した時に、雲だらけで「あーー」ってなったことありますよね。

Sentinel-2 で一番簡単に雲除去ができる方法をご紹介します。

結論

データ取得の時に SCL.jp2 も一緒に取得して以下のコマンド

gdal_calc.py -A input.tif -B R20m/SCL.jp2 --outfile=result.tif \
	--calc="A*(1 - logical_or(B==8, B==9))" --NoDataValue=0

解説

バンド演算

https://gdal.org/programs/gdal_calc.html

このコマンドによってバンドやファイル同士の演算が可能になります。

Sentinel-2 Scene Classification Map

https://custom-scripts.sentinel-hub.com/custom-scripts/sentinel-2/scene-classification/

クラスのインデックス (8, 9) は薄雲と雨雲なのでこの情報を利用します。

シーン選定

フィリピンの画像で雲を透過してみます

https://apps.sentinel-hub.com/eo-browser/?zoom=12&lat=10.32927&lng=123.8166&themeId=DEFAULT-THEME&visualizationUrl=https%3A%2F%2Fservices.sentinel-hub.com%2Fogc%2Fwms%2Fbd86bcc0-f318-402b-a145-015f85b9427e&datasetId=S2L2A&fromTime=2023-12-10T00%3A00%3A00.000Z&toTime=2023-12-10T23%3A59%3A59.999Z&layerId=1_TRUE_COLOR&demSource3D="MAPZEN"

以下で取得します。

aws s3 cp s3://sentinel-s2-l2a/tiles/51/P/WM/2023/12/10/0/ ./S2/ --recursive --request-payer requester

サンプル可視化

RGB画像ではこのようになっています。

すると、 R20m の配下に SCL.jp2 に雲の情報が含まれています。
なので同じく R20m/B02.jp2 の青色🟦波長の画像で検証してみます。

雲除去

以下のコマンドで雲除去をします。

gdal_calc.py -A R20m/B02.jp2 -B R20m/SCL.jp2 --outfile=result.tif \
	--calc="A*(1 - logical_or(B==8, B==9))" --NoDataValue=0

可視化

結果を緑色🟢で可視化します。

ちゃんとコマンド1発で青色の雲の部分が透過できていますね。

詳細の画像

精密な除去

これに関しては、各個人でデータを設定してモデルを組んで処理するが良いと思います。
実データを使った雲検知のディープラーニングモデル開発の仕方の記事もコード付きで執筆しているので参考になれば幸いです。

https://sorabatake.jp/33451/

環境構築

GDAL は環境構築に失敗することが多いと思います。
以下の Docker 利用で簡単に作成できるのでおすすめです。

https://github.com/syu-tan/sorabatake-004

おまけ

こちら以外にも記事執筆をしているのでご参考になれば幸いです

https://zenn.dev/syu_tan

衛星データ解析として、宙畑のライターもしています。

https://sorabatake.jp/?s=秀輔

SAR 解析をよくやっていますが、画像系AI、地理空間や衛星データ、点群データに関心があります。
勉強している人は好きなので楽しく絡んでくれると嬉しいです。

https://twitter.com/emmyeil

Discussion