Google Earth Engine でLandsat画像から雲の少ない自然色画像を作成する

# はじめに
本記事では、Google Earth Engine（GEE）を使って、**特定の国（土地）を対象に、雲や雲影の影響を取り除いた自然色の衛星画像を合成・出力する方法**を紹介します。

対象国はエチオピア、対象期間は2017年の1年間とし、Landsat 8の衛星画像を利用します。

本記事に含まれるクラウドマスキングや反射率スケーリングなどの基本処理も行うことで、
対象国の全土をカバーした「雲なし自然色画像ベースマップ」を作成することができます。

作成した画像は、報告書用の地図データ、土地被覆解析などの前処理素材として活用できます。

![Ethiopia_L8_2017_SR](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4032032/c261f99d-dc91-46fb-a613-b56c8f6e7260.png)


# ゴール
このプロジェクトのゴールは、以下の2点です：

- 対象国の自然色画像（R, G, B）を、雲の影響を除去した状態で合成・描画する
- 合成画像を GeoTIFF 形式でエクスポートし、Google Drive に保存する

## 処理対象の条件：
| 項目        | 内容                                                 |
| --------- | -------------------------------------------------- |
| **対象国**   | エチオピア（Ethiopia）                                    |
| **対象期間**  | 2017年（2017年1月1日〜2017年12月31日）                       |

# 分析環境と使用データセット
| ソフト／データセット                | バージョン・URL                                                                                          |
| ----------------------- | -------------------------------------------------------------------------------------------------- |
| **Google Earth Engine** | [Code Editor（ブラウザ版）](https://code.earthengine.google.com/) |
| **国境ポリゴン**              | `USDOS/LSIB_SIMPLE/2017`<br>[LSIB 2017: Large Scale International Boundary Polygons, Simplified](https://developers.google.com/earth-engine/datasets/catalog/USDOS_LSIB_SIMPLE_2017?hl=ja) |
| **Landsat 8　画像**              |  `LANDSAT/LC08/C02/T1_L2`<br>[USGS Landsat 8 Level 2, Collection 2, Tier 1](https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1_L2?hl=ja)|

# 手順

## 1. 対象領域のインポート（国境ポリゴンの抽出）

まずは解析対象とする**エチオピアの国境ポリゴン**をGoogle Earth Engine上に読み込みます。
GEEでは、あらかじめ用意された「国境データセット」から任意の国をフィルタで抽出できます。

```php:javascript
// 1. エチオピアの国境ポリゴンを取得
var ethiopia = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
                   .filter(ee.Filter.eq('country_na', 'Ethiopia'));
```
### 解説
1行目：`ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017') `は、Google Earth Engine（GEE）の公式カタログに登録されている国際境界データセット[**「LSIB 2017：Large Scale International Boundary Polygons, Simplified」**](https://developers.google.com/earth-engine/datasets/catalog/USDOS_LSIB_SIMPLE_2017?hl=ja)を読み込んでいます。

- カタログページ
[LSIB 2017 – Large Scale International Boundary Polygons, Simplified](https://developers.google.com/earth-engine/datasets/catalog/USDOS_LSIB_SIMPLE_2017?hl=ja)

![LSIB 2017](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4032032/baa8b48f-52ee-451c-a1e1-c2795bc184c1.png)

:::note info
`USDOS/LSIB_SIMPLE/2017`の「2017」は、米国務省（USDOS）が公開した**2017年版の境界データ**を指します。
**現時点（2025年）でも、GEEで配布されている最新版**として使われています。

最新の行政区画や、より詳細な境界が必要な場合（例：スーダン・南スーダンのような係争地域）は、FAO GAUL や Natural Earth などの代替データを検討します。
:::

2行目：`.filter(ee.Filter.eq('country_na', 'Ethiopia'));`では、属性フィールド ` 'country_na' `（＝国名）が ` 'Ethiopia' `のポリゴンだけを抽出しています。

`'country_na'`は、LSIB データに含まれる英語表記の国名です（例：Kenya、Japan など）。他国を抽出する場合は、国名の箇所を書き換えます：
```php:js
.filter(ee.Filter.eq('country_na', 'Kenya'));
```
:::note alert
GEEのJavascriptでは**大文字・小文字が区別される**ため、
 `'ethiopia' `、 `'japan' `（全て小文字）ではマッチせずエラーとなります。
正しくは `'Ethiopia' ` 、`'Japan' `と入力します。
:::


これで解析対象エリア（AOI）の準備が完了です。

このステップで取り込んだポリゴン（`'Ethiopia' `）は、次以降のステップで、対象国と重なるLandsat画像だけを抽出するためのフィルタに使用したり、出力時の領域の指定に使用します。

次は衛星画像の取得へ進みます。

:::note tip
AOI（Area of Interest）は「解析対象エリア」の略称です。
GEEやGIS解析における処理や出力範囲の基準となる領域を示します。
:::

## 2. Landsat 8画像の取得

次に、**Landsat 8 の画像コレクション**から、**対象期間に撮影された画像**のうち、先ほどの**ポリゴンと重なる画像だけ**を取得します。

```php:javascript
// 2. Landsat 8 Collection 2, Level-2 SR 画像（2017年）を取得
var l8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
            .filterDate('2017-01-01', '2017-12-31')
            .filterBounds(ethiopia);
```
### 解説：
1行目：`ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')`は、読み込みたいデータセットを指定しています。

GEE公式カタログに登録されている[**「Landsat 8 Level-2, Collection 2, Tier 1（Surface Reflectance）」**](https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1_L2?hl=ja)の画像コレクションを読み込むように指示しています。

- カタログページ
[Landsat 8 Level-2, Collection 2, Tier 1（Surface Reflectance）](https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1_L2?hl=ja)

![Landsat 8 Level-2, Collection 2, Tier 1](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4032032/9ac73883-dc50-4e54-902a-988bac5ba56d.png)

2行目以降では、フィルタ条件を設定し、対象画像を絞り込んでいます。

| 処理                        | 内容                                                               |
| ------------------------- | ---------------------------------------------------------------- |
| ` .filterDate('2017-01-01', '2017-12-31')`        | 2017年1月1日〜12月31日までの**1年間の画像に限定**|
| `.filterBounds(ethiopia);` | 抽出済みの**エチオピア国境ポリゴンと重なる画像だけ**に絞り込み|

:::note info
このデータセット（Level-2 SR）は**大気補正済みの反射率データ**です。
NDVIやNDWIなどの指標計算や、雲・雲影を除去する処理に最適です。
:::

## 3. 雲マスクと合成処理（雲の少ない１枚画像を作成）

Landsat画像には、雲や雲影（雲による影）が含まれており、そのままでは解析や可視化の際に支障が出る可能性があります。特に、モニタリングや土地被覆分類などでは地表の状態が見えていることが重要なため、雲の除去処理は必須です。

このステップでは、以下の3つの処理を行い、**1年間の画像から雲の少ない自然色画像を合成**します。

 3-1\. Reflectance（反射率）への変換
 3-2\. 雲マスク処理
 3-3\. 中央値合成（1年分の画像を1枚に統合）

### 3-1.Reflectance（反射率）への変換

衛星画像の値（DN値）を、反射率（Reflectance）に変換します。

```php:javascript
// 3-1. 反射率への変換（SRスケーリング）
function prepSR(img) {
  var sr = img.select('SR_B.*')               // B1〜B7の反射率バンドをまとめて取得
              .multiply(2.75e-05)             // スケール係数で割り戻し
              .add(-0.2)                      // オフセット補正
              .clamp(0, 1);                   // 反射率を0〜1の範囲に制限
```
#### 解説：
反射率（Reflectance）とは、地表が太陽光をどれくらい反射しているかを示す値で、通常は 0〜1 の実数で表されます（例：0.13 = 13%反射）。

しかし、Landsat 8 の元データ（Level-2 SR）は、各ピクセルに0〜10,000の整数（DN値）でスケーリングされた値が入っており、このままでは可視化や解析（NDVIなど）に利用できません。

そのため、上記のように 解析等で使える、意味のある実数の反射率に変換します。

:::note tip
**スケール係数、オフセットの根拠**
NASAやUSGSの公式仕様に基づき、以下のスケーリング式が定義されています：
- スケール係数：`0.0000275`
- オフセット：`-0.2`

例：DN値 = 5000の場合
　→ 5000 × 0.0000275 - 0.2 = 0.1375（＝13.75%）
上記により、意味のある反射率が算出されます。
:::

### 3-2.雲マスク処理

反射率に変換した画像には、まだ雲や雲影（＝雲によってできた影）が含まれています。
そのため、各ピクセルの状態をもとに、雲と影を自動で検出し、除去する処理（マスク処理）を行います。

```php:javascript
// 3-2. 雲マスク処理
  var qa = img.select('QA_PIXEL');  // 各ピクセルの状態を示すQAバンドを取得

  var mask = qa.bitwiseAnd(1 << 3).eq(0)  // Bit 3 = 雲影 → 0（雲影なし）
               .and(qa.bitwiseAnd(1 << 4).eq(0)); // Bit 4 = 雲 → 0（雲なし）

  return sr.updateMask(mask)
           .copyProperties(img, ['system:time_start']);
}
```
#### 解説：
雲や雲影の情報は、Landsat画像に付属している`QA_PIXEL`バンドに格納されています。これは各ピクセルの品質情報をビット形式（フラグ）で表したもので、「このピクセルは雲か？影か？」を判別するために使える補助データです。

この処理では、雲影（Bit 3）と雲（Bit 4）が立っていない（＝0）ピクセルだけを残すことで、雲の影響が少ないピクセルだけを可視化・解析に利用できるようにしています。
 
| 処理内容              | 意味                                |
| ----------------- | --------------------------------- |
| `1 << 3`          | 「8」＝Bit 3（雲影フラグ）を意味する|
| `1 << 4`          | 「16」＝Bit 4（雲フラグ）を意味する|
| `bitwiseAnd(...)` | 指定ビットが「1」かどうかをチェックする演算 |
| `.eq(0)`          | ビットが立っていない（＝雲・雲影がない）ピクセルのみ残す |
| `.updateMask()`          | マスク処理。雲や影があるピクセルを透明にして、以後の処理で無視させる|

:::note tip
`QA_PIXEL`バンドのビット情報は、衛星データごとに異なります。
Landsatでは Bit 3＝雲影、Bit 4＝雲ですが、SentinelやMODISでは別のビット位置になります。
ビット構造を確認したい場合は、GEE公式カタログ内の「Bands」セクションをチェックします。
:::

この処理によって、雲や雲影を自動的に除去した画像を得ることができます。
次に、こうして前処理された画像を1枚にまとめる「合成処理」を行います。

### 3-3.中央値合成（1年分の画像を1枚に統合）

マスク処理により、雲のない画像が1枚ずつ作成されました。
ここでは、それら複数の画像（今回は2017年の1年分＝20〜30枚程度）を1枚の画像に統合します。

```php:javascript
// 3-3. 雲除去済み画像の中央値コンポジット
var composite = l8.map(prepSR)
                  .select(['SR_B4', 'SR_B3', 'SR_B2'])  // R, G, Bの自然色バンド
                  .median()　       //中央値で合成
                  .clip(ethiopia);  //国境ポリゴンでトリミング
```
#### 解説：
composite は、1年間の雲のないピクセルを集めた「自然色画像」です。
ここでは`.median() `関数を使って、各ピクセル位置の中央値を採用しています。

また、`.clip(ethiopia)` をつけることで、合成画像を国境ポリゴンできれいに切り抜くことができます。
`.filterBounds(ethiopia)`は、「エチオピア国境に一部でも重なるすべてのLandsatタイル画像」を取得する処理です。
そのため、国境線を少しはみ出した画像も含まれてしまいます。
このような場合は、`.median() `の後に `.clip(ethiopia)`を追加することで、画像を国境ポリゴンに沿ってトリミング（切り抜き）し、視覚的にも整った国単位の画像を得ることができます。


| 処理内容              | 意味                                |
| ----------------- | --------------------------------- |
| `.map(prepSR)`          | 画像コレクション内のすべての画像に対して、前処理関数（Reflectance変換＋雲マスク）を適用|
| `.select(...)`          | 自然な色を再現するため、赤（B4）・緑（B3）・青（B2）バンドだけを使用|
| `.median()` | 有効な（マスクされていない）値の中央値をピクセルごとに採用して合成 |

:::note tip
**なぜ「中央値（median）」を使うのか？**
平均（mean）を使うと、雲（= 明るい）や影（= 暗い）などの極端な値に引きずられやすく、画像が不自然になります。
中央値（median）であれば、外れ値の影響を受けにくく、自然で安定した見た目の画像が得られます。
:::

このステップで得られた`composite`という変数には「雲と雲影を除去済み」で「複数のLandsat画像（SR）を統合」した「自然色のベースマップ」の画像ができました。

この画像は、次のステップで 画面に表示したり、GeoTIFF形式で出力して利用できます。

## 4. 表示と出力（画面への可視化とGeoTIFF書き出し）

これまでのステップで完成したcomposite画像は、雲や雲影が除去された反射率画像であり、1年間のLandsat画像を統合した自然色の1枚となっています。

このステップでは、それを Google Earth Engine（GEE）のマップ上に表示し、GeoTIFF形式でGoogle Driveへ出力できるようにします。

### 4-1. 画面への表示
```php:javascript
// 4-1. 画面表示（自然色）
Map.centerObject(ethiopia, 6);  // エチオピアを中心にズーム
Map.addLayer(composite,         // 表示する画像（R,G,B合成）
             {min: 0.05, max: 0.4},
             'Ethiopia L8 2017');
```

#### 解説：
| 処理                      | 内容                            |
| ----------------------- | ----------------------------- |
| `Map.centerObject(...)` | エチオピアのポリゴンを画面中央にズーム（ズームレベル 6） |
| `Map.addLayer(...)`     | 合成画像を自然色（RGB）としてマップに重ねて表示     |


:::note tip
`min`と`max`の値は、反射率画像の明るさを調整する表示設定です。
一般的にSRバンドの可視化には min: 0.05, max: 0.4 が自然な色調になります。
地形や土地被覆によって微調整すると見やすさが向上します。
:::

![国境ポリゴンで切り抜いた合成画像](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4032032/ab0e0f9b-b0c9-47f5-9db6-bd7e73c9df5c.png)

### 4-2. GeoTIFFとして出力（任意）
```php:javascript
// 4-2. GeoTIFF出力
Export.image.toDrive({
  image: composite.unmask(0),   // 書き出す画像
  description: 'Ethiopia_L8_2017_SR',  // ファイル名（拡張子は自動）
  region: ethiopia.geometry(),         // 出力範囲（国境ポリゴンでクリップ）
  scale: 30,                     // 解像度（30m/pixel）
  maxPixels: 1e13       　　　　   // 上限ピクセル数（大きな国土でも対応）
});
```
#### 解説：
このコードを使うことで、GeoTIFF形式の衛星画像をGoogle Driveに保存できます。
後からQGISやArcGISなどで読み込んで使うことが可能です。

書き出したGeoTIFFは、QGISやArcGISといったGIS環境で読み込んで利用できます。
また、報告書・地図作成・重ね合わせ解析などにも活用できます。

| パラメータ         | 説明                             |
| ------------- | ------------------------------ |
| `image`       | 出力対象の画像（ここでは `composite`）      |
| `description` | 出力ファイルの名前（拡張子 `.tif` は自動で付きます） |
| `region`      | 書き出す範囲（ここでは `ethiopia` の国境）    |
| `scale`       | 空間解像度（Landsat標準の 30m/pixel）    |
| `maxPixels`   | 書き出し上限ピクセル数。広域を扱うときは必ず指定       |

以上で、対象国全体の雲を除いた自然色ベースマップをGEEで作成し、表示・出力まで完了しました。

# まとめ
ここまでの処理をすべてまとめた Google Earth Engine（GEE）用の完成コードは以下のとおりです。

これをコピーして、GEEのブラウザ環境で実行すれば、エチオピア全土の雲を除去した自然色画像を簡単に生成できます。

## コード全文

```php:javascript
// 1. エチオピアの国境ポリゴンを取得
var ethiopia = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
                   .filter(ee.Filter.eq('country_na', 'Ethiopia'));
                   
// 2. Landsat 8 Collection 2, Level-2 SR 画像（2017年）を取得
var l8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
            .filterDate('2017-01-01', '2017-12-31')
            .filterBounds(ethiopia);

// 3. 雲マスクと合成処理
// 3-1. 反射率への変換（SRスケーリング）
function prepSR(img) {
  var sr = img.select('SR_B.*')               // B1〜B7の反射率バンドをまとめて取得
              .multiply(2.75e-05)             // スケール係数で割り戻し
              .add(-0.2)                      // オフセット補正
              .clamp(0, 1);                   // 反射率を0〜1の範囲に制限
// 3-2. 雲マスク処理
  var qa = img.select('QA_PIXEL');  // 各ピクセルの状態を示すQAバンドを取得

  var mask = qa.bitwiseAnd(1 << 3).eq(0)  // Bit 3 = 雲影 → 0（雲影なし）
               .and(qa.bitwiseAnd(1 << 4).eq(0)); // Bit 4 = 雲 → 0（雲なし）

  return sr.updateMask(mask)
           .copyProperties(img, ['system:time_start']);
}
// 3-3. 雲除去済み画像の中央値コンポジット
var composite = l8.map(prepSR)
                  .select(['SR_B4', 'SR_B3', 'SR_B2'])  // R, G, Bの自然色バンド
                  .median()　       //中央値で合成
                  .clip(ethiopia);  //国境ポリゴンでトリミング
                  
// 4. 表示と出力
// 4-1. 画面表示（自然色）
Map.centerObject(ethiopia, 6);  // エチオピアを中心にズーム
Map.addLayer(composite,         // 表示する画像（R,G,B合成）
             {min: 0.05, max: 0.4},
             'Ethiopia L8 2017');
// 4-2. GeoTIFF出力
Export.image.toDrive({
  image: composite.unmask(0),   // 書き出す画像
  description: 'Ethiopia_L8_2017_SR',  // ファイル名（拡張子は自動）
  region: ethiopia.geometry(),         // 出力範囲（国境ポリゴンでクリップ）
  scale: 30,                     // 解像度（30m/pixel）
  maxPixels: 1e13       　　　　   // 上限ピクセル数（大きな国土でも対応）
});
```

## 実行方法（GEE Code EditorでRunする方法）
 1\. [Google Earth Engine Code Editor](https://code.earthengine.google.com/)にアクセス
 2\.エディタ画面に上記コードをすべてコピー＆ペースト
 3\.エディタ左上の▶「Run」ボタンをクリック
　 → 画面下のマップエリアに、エチオピアの雲を除去した自然色画像が表示されます
 4\.GeoTIFFを書き出す場合は、画面右枠内の「Tasks」タブを開く
 - 「Tasks」タブのUnsubmitted tasks（未送信タスク） に今回作成したファイル名（descriptionで設定したファイル名）が表示されます
 - ファイル名横にある「Run」ボタンをクリック
 - 「Run」ボタンをクリックすると、「Task: Initiate image export」という設定ダイアログが開くので、ダイアログ内の「Run」ボタンクリック。エクスポートが開始されます
 5\.数分〜十数分後、GeoTIFFファイルが次の場所に保存されます：
   Google Drive > 「earthengine」フォルダ（自動作成） > `[ファイル名].tif  `

このコードは、年・対象国・解析指標を変えるだけで、教育や研究、さまざまな実務で応用できます。

# 参考資料

- https://qiita.com/ishimotojuri/items/cb6d572392cdb604c4c7
- https://gunmagisgeek.com/blog/googleearthengine/6238