雲のない画像を入手するには？

はじめに

先日にQGIS上で衛星画像を表示する方法を紹介しました（先日の記事はこちら）。

ただ、各プロパイダーのウェブサイトから衛星画像を入手する際に、できる限り雲のない画像をダウンロードしたいものの、検索要件に当てはまった画像のリストから一枚ずつ雲の割合を確認する必要があり、時間がかかりました。。

他方で、Google Earth Enginという無料のプラットフォームでは、衛星画像の収集が簡単に行えるため、今回はGoogle Earth Engineを使って雲の影響が少ない画像の入手・表示方法を紹介したいと思います。

分析環境

今回はGoogle Earth Engineを使いました。

手順

1. 関心領域のインポート（今回はコンゴ民主共和国とします）
2. Landsat画像の入手
   2.1. 中央値による雲なし合成画像の入手
   2.2. マスク処理による雲なし合成画像の入手

1. 関心領域のインポート

"country_co"はここから確認できます。

sample.js

var DRC = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
         .filter(ee.Filter.eq("country_co", "CG"));

Map.centerObject(DRC, 5);
Map.addLayer(DRC, {}, "DRC");

2. Landsat画像の入手

雲のない画像を探すには、大きく分けて以下の3つの方法があります。

1. 一枚の画像を雲の割合が低い順に並べて、最も雲のない画像を選ぶ
2. ある期間に撮影された複数の画像を中央値で合成する
3. 雲によって隠れてしまっている領域に雲がかかっていない衛星画像のデータを合成する（マスク処理）

ただし、上記のうち1は一枚の画像のみで示される狭い領域で使える手法であり、今回のような複数の画像で表示される広い領域では、使えません。そのため、今回は上記の2と3を使用して雲のない画像を表示させたいと思います。

2.1. 中央値による雲なし合成画像の入手

sample.js

var L8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_RT_TOA")
           .filterDate('2022-01-01','2022-12-31')
           .filterBounds(DRC)
           .median()
           .clip(DRC);

var vizTrue = {
  bands: ['B4', 'B3', 'B2'],
  min: 0,
  max: 0.4,
};

Map.addLayer(L8, vizTrue, 'L8 True Color DRC');

今回はコンゴ民主共和国を対象としているのですが、熱帯雨林地域なので一年を通して雲が多く、また対象領域が国レベルで比較的大きいため、雲のある画像が入ってしまっています。

2.2. マスク処理による雲なし合成画像の入手

次にマスク処理を行い、雲によって隠れてしまっている領域に雲がかかっていない衛星画像のデータを合成して、一枚の画像を作成します。

マスク処理は、特定の部分のみを表示（抽出）し、それ以外の部分を表示しない（黒色画像または白色画像）ようにする画像処理のことをいいます。（出所：リモートセンシング画像処理）

上記のLandsat 8のデータには、QA_PIXELというピクセルの品質評価（Quality Assessment: QA）を示すバンドが含まれています。このバンドには対象のピクセルに雲や雲の影が含まれているかを示すデータが含まれています。例えば、QA_PIXEL(雲=3->)であれば、そのピクセルは雲の情報が含まれているということを意味します。

このQA_PIXELには情報の種類によってビット（情報量の最小単位。二進法（0 or 1）の数字）が割り当てられています。0であればその情報が入っていない、1であればその情報が入っている可能性が高い、ということです。

各ビットが示す情報は以下の通りです。
ピクセル品質評価ビット指数
ビット 1: 水蒸気雲
ビット 2: 巻雲（うっすらとした白い雲）
ビット 3: 雲
ビット 4: 雲の影
ビット 5: 雪
ビット 6: 快晴（雲や巻雲がない状態）
ビット 7: 水
ビット8～9：雲の信頼度
ビット10～11: 雲の影の信頼度
ビット12～13: 雪の信頼度
ビット14～15: 巻雲の信頼度
（出所：Landsat Collection 2 Pixel Quality Assessment Bit Index）

オンオフ（0 or 1）を確認することで、画像内の雲の領域を判別しマスクを作成することができます。

まずはfunctionという新たな関数を作成します。functionの流れは以下の通りです。
①上記のビット 3: 雲、ビット 4: 雲の影の変数を作成します（そのピクセルに雲があれば1、なければ0、雲の影があれば1、なければ0）。
②各画像からQA_PIXELという情報（バンド）だけを抽出する変数を作成し、それをqaという名前にします。
③②の情報に関して、①で作成したビット3：雲が0で、かつビッド4：雲の影も0であれば1を取る変数を作成し、それをmaskという名前にします。
④③で作成したmaskの情報に基づいて画像のデータを更新します（そのピクセルのmaskが0である場合、ピクセルの情報を0に返す（その部分を表示しない＝黒色画像）、１である場合ピクセルの情報を保持（その部分を表示する）する。

sample.js

var mask_clouds_landsat8 = function(image) {
  // Bits 3 and 4 are cloud and cloud shadow, respectively.
  var cloudBitMask = (1 << 3);
  var cloudShadowBitMask = (1 << 4);
  
  // Get the pixel QA band.
  var qa = image.select('QA_PIXEL');
  
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
               .and(qa.bitwiseAnd(cloudShadowBitMask).eq(0));
  
  // Mask image with clouds and shadows
  return image.updateMask(mask);
};

また、画像をすべての領域分（世界中）表示させようとすると時間がかかるので、画像をコンゴ民主共和国のみにして、表示させる関数を作成しておきます。

sample.js

var clip_DRC = function(image) {
  return image.clip(DRC);
};

2022年1月1日から2022年12月31日までのすべての画像に対し、上記で作成した関数を適用させ、最終的にすべての画像のピクセルの中央値(median)を抽出した一枚の画像を作成します。

sample.js

var L8_masked_median = ee.ImageCollection('LANDSAT/LC08/C02/T1_RT_TOA')
                       .filterDate('2022-01-01','2022-12-31')
                       .filterBounds(DRC)
                       .map(clip_DRC)
                       .map(mask_clouds_landsat8)
                       .median();

Map.addLayer(L8_masked_median,vizTrue,'Landsat 8 - RGB - DRC Masked Median Composite');

マスク処理を行うと、中央値による合成画像よりは雲がなくなり、きれいになったように思います。

コンゴ民主共和国の首都キンシャサを拡大し、航空画像と比較したものが下図です。
上：航空写真、下：Landsat画像（マスク処理したもの）

特にLandsat画像の方が見にくく、何があるのかわかりにくいため、植生は赤色で示されるフォースカラー合成で表示してみます。

トゥルーカラー合成の場合はR（赤）：赤（Landsat-8　バンド4）、G（緑）：緑（バンド3）、B（青）：青（バンド2）に相当するバンドのデータを割り当てるのに対し、フォースカラー合成の場合は、R（赤）：近赤外（バンド5）、G（緑）：赤（バンド4）、B（青）：緑（バンド3）に相当するバンドのデータを割り当てて表示します。

※バンドの種類と色はこちらから確認できます。

sample.js

var vizFalse = {
  bands: ['B5', 'B4', 'B3'],
  min: 0,
  max: 0.5,
};

Map.addLayer(L8_masked_median,vizFalse,'Landsat 8 - False - DRC Masked Median Composite');

このあたりは、中央を除きほぼ緑に囲まれており、右上にはコンゴ川（フォースカラーでは黒色）があることがわかります。

おわりに

ただ中央値で合成画像を作成するよりは、マスク処理を行い合成画像を作成する方が相対的に雲のないきれいな画像を入手できたと思います。熱帯地帯で雲のない一枚の画像を入手することは非常に困難であり、相対的に雲のない画像を探そうとも、プロバイダーのサイトから地道に目視で確認していくのには非常に時間がかかります。

今回ご紹介したようにGoogle Earth Engineを使うことで時間の節約にもなるので、ぜひ一度試してもらえたらと思います。

また今回はLandsatの画像のみでしたが、次はSentinelの画像でマスク処理を行ってみたいと思います。