【Python】衛星データから地図上で植生状態を確認する方法【NDVI】

【Python】衛星データから地図上で植生状態を確認する方法【NDVI】

衛星画像を利用した解析を行いたいと思ったことはありませんか?衛星画像を利用したリモートセンシングは、温暖化や自然災害などの環境変動の調査、都市開発やインフラの監視といった様々な分野で活用されています。

現在、NASAのLandsatやESAのSentinelといった衛星のデータが無料で公開されています。そのデータを使えば、都市や森林といった地表面の様子が誰でも調べることができます。

今回はSentinel-2という衛星の画像データとGoogle Earth Engineというサービスを用いて、森林や田畑などの植生を示す指標であるNDVIを地図にプロットするところまでを扱います。

動作検証済み環境
macOS Monterey(12.6), python3.9.13, earthangine-api 0.1.325, colorcet 3.0.0, folium 0.12.1.post1

NDVI(正規化植生指数)について

人工衛星のひとつである地球観測衛星は、地表面に何が存在するかという情報をさまざまなセンサを利用して観測しています。
観測センサの一つである光学センサは、地表面からの反射光を撮影し、可視光線から近赤外線の各波長の反射強度の違いから、地表面に何が存在するのかを判別します。

地球の地表面に存在する植物、土壌、水(海)の反射スペクトルは左に示すグラフのような傾向をしています。
中でも、植物は特徴的なスペクトルを示しています。植物は青色光や赤色光のエネルギーで光合成をおこなっています。
そのため、植物の反射スペクトルをとると400 nm付近の青色光や600から700 nm付近の赤色光の反射強度が低くなります。
反対に700 nmより長波長の近赤外線は反射するという特性があるため、植物の生育状態を調べる植生調査に応用されています。

植生調査では、下記に示すようなNDVI(Normalized Difference Vegetation Index:正規化植生指数)という指標が利用されています。
この指標を用いることで、土壌や海水、建物などの人工物と植生の状態が判別できます。

$$
NDVI = NIR-RED/NIR+RED
$$

$NIR$:近赤外線の反射率、$RED$:赤色光の反射率

今回は無料で利用可能である点、比較的分解能が高い点、赤色光と近赤外線の光学センサを積んでいる点から、ESAのSentinel-2という衛星のデータを利用します。

Google Earth Engineについて

Google Earth Engine(以下GEE)はGoogleのクラウドやAPIを利用して衛星画像データの入手や解析が可能なサービスで、研究や教育活動では無料で利用することができます。

これまで人工衛星のデータ利用は、衛星を所有する各機関への利用登録、衛星画像データのダウンロード、専用のソフトウェアによる描画と解析と多くの手間が衛星画像1枚ごとに必要でしたが、GEEはすべてGoogleのクラウド上で完結することができます。
今回はGEEを利用するため、利用登録方法を紹介します。

GEEの利用登録まで

GEEの公式サイトにアクセスし「Sign Up」をクリックします。クリックしたらGoogleアカウントのログインが求められるので、適当なアカウントでログインします。

次に利用者の名前・所属・所属機関の区分・居住国・利用用途を入力し、利用規則を確認後、「SUBMIT」をクリックします。
GEEは研究や教育活動では無料ですが、GEEを利用したアプリ販売などの商用利用は有償利用登録が必要となるので注意が必要です。

利用登録が完了したらGoogleから「Welcome to Google Earth Engine!」という件名のメールが届くので確認後利用できます。

GEEを使ってNDVIを地図上にプロットする

GEEはJavaScriptまたはPythonでプログラミングできます。
JavaScriptはCode Editorというブラウザ上で操作可能なツールが用意されているため直感的に操作がしやすく、Pythonは機械学習への連携が容易といったメリットがあります。
今回はPythonでコードを書きます。

GEE利用の準備(Python)

pip を使って今回使用するパッケージをインストールします。

$ pip install earthengine-api  # GEEのパッケージ
$ pip install colorcet  # カラーバーを作成するパッケージ
$ pip install folium  # 地図を作成するパッケージ

Successfullyという文言が表示されたらインストール完了です。

ソースコード

import colorcet as cc
import ee
import folium

# GEEの認証
ee.Authenticate()
ee.Initialize()

# foliumマップにGEEを表示させる関数
def add_ee_layer(self, ee_image_object, vis_params, name):
    map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)
    folium.raster_layers.TileLayer(
        tiles=map_id_dict['tile_fetcher'].url_format,
        attr='Map Data © <a href="https://earthengine.google.com/">Google Earth Engine</a>',
        name=name,
        overlay=True,
        control=True
    ).add_to(self)

# 雲に覆われているピクセルの補正処理関数
def maskS2clouds(image):
    qa = image.select('QA60')
    cloudBitMask = 1 << 10
    cirrusBitMask = 1 << 11
    mask = qa.bitwiseAnd(cloudBitMask).eq(0) and qa.bitwiseAnd(cirrusBitMask).eq(0)
    return image.updateMask(mask).divide(10000)

# NDVIの計算用関数
def calc_ndvi(image):
    return ee.Image(image.expression(
    '(NIR-RED)/(NIR+RED)', {
      'RED': image.select('B4'),
      'NIR': image.select('B8')
}))

# Sentinel-2のデータの収集
Sentinel2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')

# 冬と夏の情報を抽出
winter = Sentinel2.filterDate('2020-12-01','2021-02-28').filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20)).map(maskS2clouds)
summer = Sentinel2.filterDate('2021-07-01','2021-09-30').filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20)).map(maskS2clouds)

# NDVIと中央値の計算
ndvi_winter = winter.map(calc_ndvi)
ndvi_winter = ndvi_winter.median()
ndvi_summer = summer.map(calc_ndvi)
ndvi_summer = ndvi_summer.median()

# NDVIの夏と冬の値の差を計算
ndvi = ndvi_summer.subtract(ndvi_winter)

# マップの中心に表示する緯度経度
lat, lon = 36.0, 140.0

# カラーバーの設定
visualization = {
  'min': 0.0,
  'max': 1.0,
  'palette': cc.rainbow
}

# マップを表示した際の初期位置と拡大値の設定
my_map = folium.Map(location=[lat, lon], zoom_start=10)

# foliumマップにadd_ee_layer関数を追加
folium.Map.add_ee_layer = add_ee_layer

# foliumマップに計算したNDVIを表示させる
my_map.add_ee_layer(ndvi_winter, visualization, 'NDVI_2020winter')
my_map.add_ee_layer(ndvi_summer, visualization, 'NDVI_2021summer')
my_map.add_ee_layer(ndvi, visualization, 'NDVI_W-S')

# RGBの表示設定
visualization_RGB = {
    'min': 0.0,
    'max': 0.3,
    'bands': ['B4', 'B3', 'B2']
}

# foliumマップにRGBを表示させる
my_map.add_ee_layer(winter.median(), visualization_RGB, 'RGB_2020winter')
my_map.add_ee_layer(summer.median(), visualization_RGB, 'RGB_2021summer')

# 複数のマップを重ね合わせる
my_map.add_child(folium.LayerControl(collapsed = False).add_to(my_map))

# マップの保存
my_map.save('ndvi_test.html')

上記のコードをndvi_test.pyという名前でDesktop/LabCode/python/satellite-analysisディレクトリに保存します。

プログラムの実行

ターミナルを開き

$ cd Desktop/labcode/python/satellite-analysis

と入力し、ディレクトリを移動します。あとは以下のコマンドを実行します。

$ python ndvi_test.py

すぐにndvi_test.htmlというファイルが作成されていると思います。

実行結果

上記のプログラムで作成したndvi_test.htmlファイルを開くと、関東平野を中心とした次のような地図が表示されるはずです。
この状態では右上のレイヤーのチェックボックスがすべてチェックされているため、一番下のRGB_2021summerが表示された状態です。

次にチェックボックスのNDVI_2021summerのみにチェックをいれた状態にすると、次のような地図が表示されます。
今回作成したカラーバーでは青→緑→黄色→赤の順にNDVIの値が高くなっています。
表示した範囲ではおそらく黄色やオレンジ色の場所が植物が繁茂している場所と判断できそうです。

最後にNDVI_W-Sの確認をします。この地図ではNDVIが冬から夏にかけて増加した場所が強調されます。
そのため落葉樹や田畑といった季節によって大きく状態が変化する植物が存在する場所が確認できます。
表示した範囲では赤の円で囲った部分に赤く強調されたところが存在します。

拡大すると次のような地図が表示できます。
この場所は印旛沼という千葉の沼で、ヒシという水面に浮遊する植物が夏場に繁茂するようです。
周りの田畑に比べて赤く表示されているのは、土壌より水のほうが赤外線を多く吸収するため、沼である印旛沼のNDVIの季節差が大きくなった影響と予想できます。

コードの解説

import colorcet as cc
import ee
import folium

ee.Authenticate() # Google認証トークンの取得と設定
ee.Initialize() # APIセッションの初期化

ここでは、今回使用するモジュールのインポートとGEEの操作をするための認証処理、APIセッションの初期化を行います。

初めてGEEの操作を行う場合、eeモジュールのAuthenticate関数を用いてGoogleの認証トークンを取得する必要があります。プログラム実行時にはGoogleへのログインとトークン発行を行います。2回目の実行以降は不要です。反対にInitialize関数は毎回実行する必要があります。

def add_ee_layer(self, ee_image_object, vis_params, name):
    map_id_dict = ee.Image(ee_image_object).getMapId(vis_params) # 表示させるGEEの画像
    folium.raster_layers.TileLayer(
        tiles=map_id_dict['tile_fetcher'].url_format, # 表示するマップタイルの設定
        attr='Map Data © <a href="https://earthengine.google.com/">Google Earth Engine</a>', # 地図に重ねる衛星データの権利者の表示設定
        name=name, # 地図に重ね合わせるデータの名前
        overlay=True, # マップタイルにレイヤーとして重ね合わせる設定
        control=True # レイヤー操作を可能にする設定
    ).add_to(self)

foliumのマップにGEEの画像データを読み込ませて表示させる関数を定義します。5行目のattrはマップデータに表示する権利者の設定です。今回は衛星画像データ入手先のGEEを設定します。7行目のoverlayと8行目のcontrolは作成したNDVIのデータを地図上でレイヤーとして表示させる設定になります。

この設定で下記のように表示させたいデータをチェックボックスとして選択できるようになります。

def maskS2clouds(image):
    qa = image.select('QA60') # 雲データマスク処理用のバンドを抽出
    cloudBitMask = 1 << 10 # 雲の情報
    cirrusBitMask = 1 << 11 # 巻雲の情報
    mask = qa.bitwiseAnd(cloudBitMask).eq(0) and qa.bitwiseAnd(cirrusBitMask).eq(0) # マスク処理用の変数の定義
    return image.updateMask(mask).divide(10000) # マスク処理した画像データを返す

雲ピクセルを補正する関数を定義します。今回利用するSentinel-2のデータには雲補正用のデータが「QA60」として用意されています。
QA60では雲の情報がBit 10とBit 11として割り当てられており、それぞれが1を取るとき雲が存在すると判断されます。
5行目ではそれぞれが0を取り雲が存在しない場合を1、それ以外の場合を0とする変数を定義します。
6行目のreturn分で雲が存在するピクセルをマスクした画像データを返します。

def calc_ndvi(image):
    return ee.Image(image.expression(
    '(NIR-RED)/(NIR+RED)', { # NDVIの計算
      'RED': image.select('B4'), # 赤色光データの抽出
      'NIR': image.select('B8') # 近赤外線データの抽出
}))

NDVI計算用の関数を定義します。Sentinel-2ではB4が665 nmの赤色光、B8が833 nmの近赤外線の情報を持っています。

# Sentinel-2のデータの収集
Sentinel2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')

# 冬と夏の情報を抽出
winter = Sentinel2.filterDate('2020-12-01','2021-02-28').filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20)).map(maskS2clouds)
summer = Sentinel2.filterDate('2021-07-01','2021-09-30').filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20)).map(maskS2clouds)

2行目のee.ImageCollection(‘COPERNICUS/S2_SR_HARMONIZED’)でSentinel-2のデータを収集します。
5行目、6行目で収集したデータから目的に合致したデータを抽出します。

今回は季節の植生変化をとらえるために冬と夏の期間を.filterDateで指定します。さらに.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))で被雲率が20%以下のデータを抽出し、.map(maskS2clouds)で雲ピクセルの補正を行います。

# NDVIと中央値の計算
ndvi_winter = winter.map(calc_ndvi)
ndvi_winter = ndvi_winter.median()
ndvi_summer = summer.map(calc_ndvi)
ndvi_summer = ndvi_summer.median()

# NDVIの夏と冬の値の差を計算
ndvi = ndvi_summer.subtract(ndvi_winter)

先ほどデータ抽出をおこなったwinterとsummerのNDVIを計算し、それぞれの中央値を計算します。
8行目では冬と夏のNDVIの中央値の差を計算します。これにより、冬から夏にかけて植生が大きく変化した場所が分かります。

# マップの中心に表示する緯度経度
lat, lon = 36.0, 140.0

# カラーバーの設定
visualization = {
  'min': 0.0,
  'max': 1.0,
  'palette': cc.rainbow
}

# マップを表示した際の初期位置と拡大値の設定
my_map = folium.Map(location=[lat, lon], zoom_start=10)

# foliumマップにadd_ee_layer関数を追加
folium.Map.add_ee_layer = add_ee_layer

# foliumマップに計算したNDVIを表示させる
my_map.add_ee_layer(ndvi_winter, visualization, 'NDVI_2020winter')
my_map.add_ee_layer(ndvi_summer, visualization, 'NDVI_2021summer')
my_map.add_ee_layer(ndvi, visualization, 'NDVI_W-S')

# RGBの表示設定
visualization_RGB = {
    'min': 0.0,
    'max': 0.3,
    'bands': ['B4', 'B3', 'B2'] # RGB作成用の赤・緑・青光データの抽出
}

# foliumマップにRGBを表示させる
my_map.add_ee_layer(winter.median(), visualization_RGB, 'RGB_2020winter')
my_map.add_ee_layer(summer.median(), visualization_RGB, 'RGB_2021summer')

# 複数のマップを重ね合わせる
my_map.add_child(folium.LayerControl(collapsed = False).add_to(my_map))

# マップの保存
my_map.save('ndvi_test.html')

最後に計算結果をマップにプロットし保存します。

最後に

今回は衛星画像データを利用した代表的な指標であるNDVIの算出と地図へのマッピングを行いました。NDVIは植生状態が確認できるので間接的に森林減少や砂漠化などがモニターできます。今回作成したHTMLファイルは場所を絞らずに全世界にプロットしているため、海外の植生状態を見てみるのも面白いかもしれません。

次回は、測定地点を絞って過去数年間のNDVI変化をグラフ化してみたいと思います。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です