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

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

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')

$ cd Desktop/labcode/python/satellite-analysis

$ python ndvi_test.py

import colorcet as cc
import ee
import folium

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

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)

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) # マスク処理した画像データを返す

def calc_ndvi(image):
    return ee.Image(image.expression(
    '(NIR-RED)/(NIR+RED)', { # NDVIの計算
      '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'] # 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')