【続編】Google Earth EngineとGoogle Colabによる人工衛星画像解析 〜無料で始める衛星画像解析（実践編）〜

はじめに

• Google Earth EngineとGoogle Colabによる人工衛星画像解析　〜無料で始める衛星画像解析（入門編）〜では、Google Earth Engine (GEE)とGoogle Colabを利用した衛星画像の取得方法を紹介しました
• 本記事では、GEEで利用できる代表的な**人工衛星データ（光学画像・植生指数・地表面温度）**の解析方法を紹介します
• 他にも様々な衛星データセットがGEEのデータカタログには保存されており、下記コードの変数（衛星名やバンド名）を変更すれば、同様の解析が行えます

衛星データを取得するコード

• 下記のコードを利用します（コードの詳細解説は前回記事をご覧ください）
• 「衛星名」「バンド名」「データの期間」「対象エリア」「保存するフォルダ名」が変数となっており、関数を実行すると指定したフォルダ（Google Drive）に衛星画像が保存されます

関数定義

# Earth Engine Python APIのインポート
import ee

# GEEの認証・初期化
ee.Authenticate()
ee.Initialize()

# GEEのデータロード
def load_data(snippet, from_date, to_date, geometry, band):
    # パラメータの条件にしたがってデータを抽出
    dataset = ee.ImageCollection(snippet).filter(
    ee.Filter.date(from_date, to_date)).filter(
    ee.Filter.geometry(geometry)).select(band)
    # リスト型へ変換
    data_list = dataset.toList(dataset.size().getInfo())
    # 対象データ数とデータリストを出力
    return dataset.size().getInfo(), data_list

# 衛星画像をGoogle Driveへ保存
def save_on_gdrive(image, geometry, dir_name, file_name, scale):
    task = ee.batch.Export.image.toDrive(**{
        'image': image,# ロードする衛星情報
        'description': file_name,# 保存するファイル名
        'folder':dir_name,# 保存先のフォルダ名
        'scale': scale,# 解像度
        'region': geometry.getInfo()['coordinates'],# 対象エリア
        'crs': 'EPSG:4326'
    })
    # Run exporting
    task.start()
    print('Done.')

変数（パラメータ）設定

## パラメーターの指定
# 衛星を指定
snippet = 'NOAA/DMSP-OLS/NIGHTTIME_LIGHTS'
# バンド名を指定
band = 'avg_vis'
# 期間を指定
from_date='2010-01-01'
to_date='2012-12-31'
# エリアを指定(日本エリアを緯度・経度を指定)
geometry = ee.Geometry.Rectangle([128.60, 29.97, 148.43, 46.12])
# 保存するフォルダ名
dir_name = 'GEE_download'
# ファイル名
file_name = 'file_name'
# 解像度
scale = 1000

処理の実行

## 処理の実行 ----------------------------------------------
num_data, data_list = load_data(snippet=snippet, from_date=from_date, to_date=to_date, geometry=geometry, band=band)
print('#Datasets; ', num_data)

## 全件保存（ファイル名は衛星IDを利用）
for i in range(data_list.size().getInfo()):
    image = ee.Image(data_list.get(i))
    save_on_gdrive(image, geometry, dir_name, image.getInfo()['id'].replace('/', '_'), scale)

##衛星データ解析① ~光学画像~

概要

• 可視光バンド（RGB）を組み合わせて表現されるカラー画像
• 一般的に利用されている「カメラ」の写真のイメージ
• RGBのそれぞれのセンサーデータをGEEから取得して、RGB合成することで可視化できる

データセット

• 衛星 = "Sentinel 2"
• バンド = TCI_R, TCI_G, TCI_B
• 解像度: 10m
• エリア: 東京（皇居を中心）

コード

• パラメータを変更して、TCI_R（可視光の赤バンド）を取得してみます

## パラメーターの指定
# 衛星を指定
snippet = 'COPERNICUS/S2_SR'
# バンド名を指定
band = 'TCI_R'

# 期間を指定
from_date='2020-04-01'
to_date='2020-04-15'
# エリアを指定(日本エリアを緯度・経度を指定)
geometry = ee.Geometry.Rectangle([139.686, 35.655, 139.796, 35.719])

# 保存するフォルダ名
dir_name = 'GEE_Sentinel_Red'
# ファイル名
file_name = 'file_name'
# 解像度
scale = 10

## 処理の実行 ----------------------------------------------
num_data, data_list = load_data(snippet=snippet, from_date=from_date, to_date=to_date, geometry=geometry, band=band)
print('#Datasets; ', num_data)

## 全件保存（ファイル名は衛星IDを利用）
for i in range(data_list.size().getInfo()):
    image = ee.Image(data_list.get(i))
    save_on_gdrive(image, geometry, dir_name, image.getInfo()['id'].replace('/', '_'), scale)

-　次にダウンロードした衛星画像を可視化してみます

# パッケージのインストール&インポート
!pip install rasterio
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import rasterio
import glob

# データの読み込み
with rasterio.open('/content/drive/My Drive/GEE_Sentinel_Red/COPERNICUS_S2_SR_20200402T012651_20200402T012651_T54SUE.tif') as src:
    arr = src.read()

# 可視化
plt.imshow(arr[0], cmap='Reds')

東京（皇居周辺）の光学画像（赤バンド）

• 中心に皇居らしきエリア、西側に代々木公園や新宿御苑が確認できます
• 次に他の緑・青バンドも取得してみます

# 緑バンドのデータ取得
## パラメーターの指定
# バンド名を指定
band = 'TCI_G'
# 保存するフォルダ名
dir_name = 'GEE_Sentinel_Green'

## 処理の実行 ----------------------------------------------
num_data, data_list = load_data(snippet=snippet, from_date=from_date, to_date=to_date, geometry=geometry, band=band)
print('#Datasets; ', num_data)

## 全件保存（ファイル名は衛星IDを利用）
for i in range(data_list.size().getInfo()):
    image = ee.Image(data_list.get(i))
    save_on_gdrive(image, geometry, dir_name, image.getInfo()['id'].replace('/', '_'), scale)

# 青バンドのデータ取得
## パラメーターの指定
# バンド名を指定
band = 'TCI_B'
# 保存するフォルダ名
dir_name = 'GEE_Sentinel_Blue'

## 処理の実行 ----------------------------------------------
num_data, data_list = load_data(snippet=snippet, from_date=from_date, to_date=to_date, geometry=geometry, band=band)
print('#Datasets; ', num_data)

## 全件保存（ファイル名は衛星IDを利用）
for i in range(data_list.size().getInfo()):
    image = ee.Image(data_list.get(i))
    save_on_gdrive(image, geometry, dir_name, image.getInfo()['id'].replace('/', '_'), scale)

• 赤バンドと同様に緑バンド、青バンドも取得することができました
• 同様にデータを読み込み可視化してみます

# データの読み込み
# Red
with rasterio.open('/content/drive/My Drive/GEE_Sentinel_Red/COPERNICUS_S2_SR_20200402T012651_20200402T012651_T54SUE.tif') as src:
    red = src.read()
# Green
with rasterio.open('/content/drive/My Drive/GEE_Sentinel_Green/COPERNICUS_S2_SR_20200402T012651_20200402T012651_T54SUE.tif') as src:
    green = src.read()
# Blue
with rasterio.open('/content/drive/My Drive/GEE_Sentinel_Blue/COPERNICUS_S2_SR_20200402T012651_20200402T012651_T54SUE.tif') as src:
    blue = src.read()

#　可視化
plt.figure(figsize=(15, 10))
plt.subplot(131); plt.imshow(red[0], cmap='Reds'); plt.title('Red')
plt.subplot(132); plt.imshow(green[0], cmap='Greens'); plt.title('Green')
plt.subplot(133); plt.imshow(blue[0], cmap='Blues'); plt.title('Blue')

東京（皇居周辺）の光学画像（赤バンド・緑バンド・青バンド）

• 最後にこの3色の画像をRGB合成して、カラー画像に変換します
• RGB合成はnumpyのdstackで結合するだけです

# RGB合成
## np.dstackで連結（Red, Green, Blueの順番に注意）
rgb = np.dstack((red[0], green[0], blue[0]))

# RGB合成した画像の可視化
plt.imshow(rgb); plt.title('RGB Image')

東京（皇居周辺）の光学画像（RGB合成）

• RGB合成することで、衛星画像らしくなりました
• 今回は特定のエリア・タイミングだけで可視化を行いましたが、人工衛星が運用されている期間内であれば任意のタイミングでデータ取得が可能です
• 季節を変えてみたり、災害前後で比較したりすると色々なものが見えてきそうです（光学画像は雲の影響を受けやすいので、雨季や天気の悪い日は画像が不明瞭になることもあります）

##衛星データ解析② ~植生指数~

概要

• 植物の分布や活性度を表す指標
• 植物に反応しやすい波長のセンサーを利用して計測
• 様々な植生指数が提案されているがNDVI（Normalized Difference Vegetation Index）が有名
• 砂漠化や森林火災の監視、都市利用分類などに利用されている

データセット

• 衛星 = "Terra(MODIS)"
• バンド = NDVI
• • 通常、NDVIは複数バンドを組み合わせた計算が必要なのですが、GEEでは計算済みのデータが保存されています（便利...!!）
• 解像度: 250m
• エリア: 東京（皇居を中心）

コード

• 各種パラメータを変更してNDVI（植生指数）を取得してみます

## パラメーターの指定
# 衛星を指定
snippet = 'MODIS/006/MOD13Q1'
# バンド名を指定
band = 'NDVI'

# 期間を指定
from_date='2005-01-01'
to_date='2005-12-31'
# エリアを指定(日本エリアを緯度・経度を指定)
geometry = ee.Geometry.Rectangle([139.686, 35.655, 139.796, 35.719])

# 保存するフォルダ名
dir_name = 'GEE_NDVI'
# 解像度
scale = 250

## 処理の実行 ----------------------------------------------
num_data, data_list = load_data(snippet=snippet, from_date=from_date, to_date=to_date, geometry=geometry, band=band)
print('#Datasets; ', num_data)

## 全件保存（ファイル名は衛星IDを利用）
for i in range(data_list.size().getInfo()):
    image = ee.Image(data_list.get(i))
    save_on_gdrive(image, geometry, dir_name, image.getInfo()['id'].replace('/', '_'), scale)

```python
## パラメーターの指定
# 衛星を指定
snippet = 'MODIS/006/MOD13Q1'
# バンド名を指定
band = 'NDVI'
# 期間を指定
from_date='2005-01-01'
to_date='2005-12-31'
# エリアを指定(日本エリアを緯度・経度を指定)
geometry = ee.Geometry.Rectangle([139.686, 35.655, 139.796, 35.719])

# 保存するフォルダ名
dir_name = 'GEE_NDVI'
# 解像度
scale = 250

## 処理の実行 ----------------------------------------------
num_data, data_list = load_data(snippet=snippet, from_date=from_date, to_date=to_date, geometry=geometry, band=band)
print('#Datasets; ', num_data)

## 全件保存（ファイル名は衛星IDを利用）
for i in range(data_list.size().getInfo()):
    image = ee.Image(data_list.get(i))
    save_on_gdrive(image, geometry, dir_name, image.getInfo()['id'].replace('/', '_'), scale)

• 同様に取得したデータを読み込み、可視化してみます

# データの読み込み
with rasterio.open('/content/drive/My Drive/GEE_NDVI/MODIS_006_MOD13Q1_2005_01_01.tif') as src:
    arr = src.read()

# 可視化
plt.figure(figsize=(20, 8))
plt.subplot(121); plt.imshow(rgb); plt.title('RGB Image'); plt.title('Optical Image')
plt.subplot(122); plt.imshow(arr[0], cmap='YlGn'); plt.title('NDVI')

東京（皇居周辺）の光学RGB画像と植生指数

• 皇居や都内の公園は植生指数が高い値になっています
• 次に取得した期間（2005年１年間）のデータを全て可視化してみます
• 一般に植生は初夏~秋にかけて活性化し、冬には衰退するので、その様子を確認してみます

# 時系列で可視化
files = glob.glob('/content/drive/My Drive/GEE_NDVI/*tif')
files.sort()

plt.figure(figsize=(15, 10))
j=0
for i in range(len(files)):
  # 画像を1シーンずつ取得して可視化
  with rasterio.open(files[i]) as src:
      arr = src.read()
  j+=1# 画像のプロット位置をシフトさせ配置
  plt.subplot(5,6,j)
  plt.imshow(arr[0], cmap='YlGn', vmin=-2000, vmax=10000)
  plt.title(files[i][-14:])# ファイル名から日付を取得
  plt.tight_layout()

東京（皇居周辺）の植生指数の変化（1月~１２月）

• 5~9月頃は植生指数が活性化している様子がみてとれます
• 若干比較がしずらいので、植生指数の合計値を算出して時系列トレンドをみてみます

# 当該エリアの　NDVIの数値を合計値を取得
sum_NDVI = []
label_NDVI = []
for i in range(len(files)):
  # 画像を1シーンずつ取得して可視化
  with rasterio.open(files[i]) as src:
      arr = src.read()
  sum_NDVI.append(arr.sum())
  label_NDVI.append(files[i][-14:-4])

# 可視化
fig,ax = plt.subplots(figsize=(15,3))
plt.plot(sum_NDVI, marker='o')
ax.set_xticks(np.arange(0,23))
ax.set_xticklabels(label_NDVI, rotation=90)
plt.show()

植生指数の変化（1月~１２月）

• 夏にかけて植生指数がピークを迎えている様子が確認できます
• 6月が大幅に現象しているのは雨季（雲）の影響で衛星が正しくデータ取得できない原因が考えられます
• 光学系の衛星は雲の影響を受けやすのが弱点ですが多様な情報を取得できる点に強みがあります（一方で雲の影響を受けないSARというタイプの衛星も存在します）
• 単年だと季節性の変化が分かりづらいので15年分のデータを取得・可視化してみました
• また移動平均も重ねて表示すると、植生指数の季節性の変化がより顕著に確認できます

植生指数の変化（15年分）

##衛星データ解析③ ~地表面温度~

概要

• "気温"が空気中の温度であるのに対して、"地表面温度"は地球の表面の温度を表す
• Land Surface Temperature(LST)と呼ばれる
• 人工構造物（建物や道路）は森林や土と比べて地表面度温度が高くなる傾向にある
• ヒートアイランド現象の監視等に利用されている

データセット

• 衛星 = "Terra(MODIS)"
• バンド = LST (Land Surface Temerature)
• 解像度: 1000m
• エリア: 東京（皇居を中心）

コード

• 各種パラメータを変更してLSTを取得してみます
• 同様に取得したデータを可視化します

## パラメーターの指定
# 衛星を指定
snippet = 'MODIS/006/MOD11A2'
# バンド名を指定
band = 'LST_Day_1km'

# 期間を指定
from_date='2005-01-01'
to_date='2005-12-31'
# エリアを指定(日本エリアを緯度・経度を指定)
geometry = ee.Geometry.Rectangle([139.686, 35.655, 139.796, 35.719])

# 保存するフォルダ名
dir_name = 'GEE_LST'
# 解像度
scale = 1000

## 処理の実行 ----------------------------------------------
num_data, data_list = load_data(snippet=snippet, from_date=from_date, to_date=to_date, geometry=geometry, band=band)
print('#Datasets; ', num_data)

## 全件保存（ファイル名は衛星IDを利用）
for i in range(data_list.size().getInfo()):
    image = ee.Image(data_list.get(i))
    save_on_gdrive(image, geometry, dir_name, image.getInfo()['id'].replace('/', '_'), scale)

# データの読み込み
with rasterio.open('/content/drive/My Drive/GEE_LST/MODIS_006_MOD11A2_2005_08_05.tif') as src:
    arr = src.read()

# 可視化
plt.figure(figsize=(20, 8))
plt.subplot(121); plt.imshow(rgb); plt.title('RGB Image'); plt.title('Optical Image')
plt.subplot(122); plt.imshow(arr[0], cmap='inferno'); plt.title('Land Surface Temperature')

東京（皇居周辺）の光学RGB画像と地表面温度

• 他データと比較して解像度が低いので若干分かりづらいです
• 該当期間を全て表示してみると、西側の方が温度が高くなる傾向にありそうです
• また地表面温度は時間解像度が高い（高頻度でデータ取得している）反面、雲の影響により正しくデータが取得できていない画像が多い点も特徴の一つです

東京（皇居周辺）の地表面温度の変化（1月~12月）

東京（皇居周辺）の地表面温度の変化（15年分）

• 植生指数と同様に長期間のデータも取得して、移動平均も含めて可視化をしてみました
• 植生と同様に季節性がありそうです（夏は高くなり、冬は低くなる）
• 解像度が低いので、もう少し広域で地表面温度を算出して、大局的なトレンドを分析した方が良さそうです

さいごに

• Google Earth EngineとGoogle Colabを用いて、代表的な衛星画像解析を紹介しました
• 上記の通りGEEとColabを利用することで、衛星名やバンド名の変数名を変更するだけで様々な衛星データの解析を行えることが分かるかと思います
• また衛星データ取得後、同一環境でPythonの便利な解析ライブラリを利用できる点が非常に便利です
• 実際の解析では、要件に応じて細かな前処理を追加したり、バイアスを除去するためにモデルを作ったりしますが、こういったサービスを利用するこで初期導入のステップは大幅に効率化できるかと思います
• また、手元の業務データをGoogle Driveへアップロードし、Colab上で衛星データとマージした解析も簡単に行えそうです
• 本記事をキッカケに様々な領域での人工衛星データ活用の一助になればと思います