Python Google Earth Engine（GEE）を用いてNDVI（植生指数）を図示する

　本記事では、Google Earth Engine（GEE）を用いて、下図のような衛星画像と 正規化植生指標（NDVI, Normalized Difference Vegetation Index）を取得する雛形コードを載せました。
　下図は鳥取砂丘の周辺の地区について、2023年の4月,8月,12月の衛星画像とNDVIの画像を図示したものです。JupyterLabを用いてインタラクティブに実行してゆきます。

▼ライブラリのインストールと、基本設定は下記リンク先の1, 2を参照ください。

hk29.hatenablog.jp

■本プログラム

はじめに、JupyterLab起動時にGEEのイニシャライズをします。

import ee
ee.Authenticate()
ee.Initialize()

次に衛星画像を取得します。抽出したい領域の緯度と経度を指定する必要があります。これは、例えばgoogleマップから取得できます。また、データ抽出期間と雲画像の除去フィルタを設定します。

import glob
import rasterio
from rasterio.plot import show
import matplotlib.pyplot as plt
from PIL import Image

################################### 入力パラメータ
# 保存するファイル名とフォルダ（グーグルドライブ内）
save_file_name = 'Tottori_Sand_Dunes'
save_folder_name = 'save_dir_for_satellite_image'

# データ抽出したい場所の矩形座標（鳥取砂丘の周辺）
coords = [[
[134.2307486167095, 35.5498136559922],
[134.2484177611878, 35.5498136559922],
[134.2484177611878, 35.5379886904566],
[134.2307486167095, 35.5379886904566],
[134.2307486167095, 35.5498136559922],
]]

# データの抽出期間
START ='2023-01-01'
END = '2023-12-09'

# 雲画像の除去フィルタ値
CLOUD_COVER_RATE = 20
###################################

# ImageCollectionの設定
# 指定できる人口衛星の種類　https://developers.google.com/earth-engine/datasets/
dataset = (
  ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED') # 人工衛星
      .filterDate(ee.Date(START), ee.Date(END)) # データ抽出の期間
      .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', CLOUD_COVER_RATE)) # 上限被覆率（雲が多い画像を除去）
      .filterBounds(ee.Geometry.Polygon(coords)) # 抽出する領域の指定 https://developers.google.com/earth-engine/guides/geometries
)

# ImageIDの取得
asset_id_list = dataset.aggregate_array('system:id').getInfo()

# 取得できたファイル数とIDを表示
print(len(asset_id_list))
asset_id_list

次に、上記ファイル名のリストから見やすいように日付を抽出します。

# データファイル名から日時を取得
date_list = []
for asset_id in asset_id_list:
    parts = asset_id.split('/')
    last_part = parts[-1]

    extracted_date = last_part[:8]
    date_list.append(extracted_date)
date_list

上記リストから、春（4月）、夏（8月）、冬（12月）のデータを抽出します。

# 必要なデータを新リストへ抽出
indices_to_extract = [4, 9, 11]
selected_elements = [asset_id_list[i] for i in indices_to_extract]
selected_elements

次に、座標参照系（CRS:Coordinate Reference System）を取得した画像ファイルから取得します。

# ImageCollectionに格納されている最初の画像の最初のバンドのCRSを取得
first_image = dataset.first()
first_band = first_image.select(0)  # 0は最初のバンドを指定
my_crs = first_band.projection().getInfo()['crs']
print('CRS:', my_crs)

下記コードでは、衛星画像を取得します。画像ファイルはgoogleドライブに保存されます。

# 衛星画像を取得する場合
my_scale = 5
my_max_pixels = 1e13

# 抽出されたファイルはGoogleドライブに保存されます
for i, asset_id in enumerate(selected_elements, start=1):
    task = ee.batch.Export.image.toDrive(
        image = ee.Image(asset_id).select(['TCI_R','TCI_G','TCI_B']), # 画像のバンド名一覧 https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR_HARMONIZED#bands
        region = ee.Geometry.Polygon(coords), # 地図領域
        folder = save_folder_name, # 保存するグーグルドライブ内のフォルダ名（もし無ければ、フォルダを作成する）
        description = f"{save_file_name}_{i}", # 保存するファイル名
        scale = my_scale, # スケールの設定（例：10メートルの解像度を示す）
        maxPixels = my_max_pixels, # 最大ピクセル数の設定
    )
    task.start()

上記を実行すると、下図のようにグーグルドライブに画像データが保存されます。

下記コードでは、NDVI画像を取得します。

# NDVI画像を取得する場合
# 抽出されたファイルはGoogleドライブに保存されます
for i, asset_id in enumerate(selected_elements, start=1):
    task = ee.batch.Export.image.toDrive(
        image = ee.Image(asset_id).normalizedDifference(['B8','B4']).rename('ndvi'), # NDVI画像を取得：B8（赤）,B4（近赤外）
        region = ee.Geometry.Polygon(coords), # 地図領域
        folder = save_folder_name, # 保存するグーグルドライブ内のフォルダ名（もし無ければ、フォルダを作成する）
        description = f"{save_file_name}_NDVI_{i}", # 保存するファイル名
        scale = my_scale, # スケールの設定（例：10メートルの解像度を示す）
        maxPixels = my_max_pixels, # 最大ピクセル数の設定
        crs = my_crs, # 測地系:座標参照系（Coordinate Reference System, CRS）
    )
    task.start()

これも、下図のようにgoogleドライブに保存されます。

※以降のコードでは、上記でgoogleドライブに保存した画像データファイルをローカルへダウンロードしてから処理します。

# NDVI画像ファイルのパスをリストで取得
files = glob.glob(f"tottori_sand/*.tif")
files = files[:3]
files

指定した衛星画像を表示します。

# 画像をプロット
for file in files:
    # 画像を開く
    image = Image.open(file)

    # 画像をMatplotlibで表示
    plt.imshow(image)
    plt.axis('off')  # 軸を表示しない設定
    plt.show()

# NDVI画像ファイルのパスをリストで取得
files = glob.glob('tottori_sand/*NDVI*.tif')
files

最後に、指定したNDVI画像を表示します。

# 画像をプロット
for file in files:
    print(file)
    ndvi_dataset = rasterio.open(file)
    ndvi_array = ndvi_dataset.read(1)
    cmap = plt.cm.RdYlBu_r # RdYlBu_r, viridis_r
    plt.imshow(ndvi_array, cmap=cmap, vmin=-1, vmax=1)
    plt.colorbar()
    plt.show()

下図のように草木がある場所は紺色、鳥取砂丘や町や道は白色、海は赤色になっている様子。

以上

<広告>

リンク

リンク