Automating GIS Processes 2024 写経 Exercise 6 （Problem 2）（１）

Problem 2: Calculating NDVI (12 points)

In this problem, we will work with the same satellite imagery data as we used during the lesson. Specifically, we will explore the Sentinel-2 satellite image from Nuuksio National Park in Finland. The data is obtained from the Copernicus Data Space Ecosystem. Our objective is to explore, process the data, and ultimately calculate the NDVI (Normalized Difference Vegetation Index).

この問題では、レッスン7と衛星画像の画像データを使います。特に、フィンランドのNuuksio国立公園のSentinel-2衛星画像を探究します。
データは Copernicus Data Space Ecosystemから取得できます。私たちの目的はデータを探究、処理、そして、究極的にはNDVI（正規化植生指標）を計算することです。

What is NDVI?

NDVI (Normalized Difference Vegetation Index) is a simple, widely used index for assessing vegetation health and coverage based on satellite imagery. NDVI is calculated by comparing the Near Infrared (NIR) and Red bands of satellite data. Vegetation reflects more NIR light and absorbs more visible red light, so the NDVI value provides insight into the density and health of vegetation.

NDVI（正規化植生指標）は単純で、衛星画像に基づき、植生の健康状態と網羅性を評価するために広く使われている指標です。
NDVIは、衛星画像のNIR（近赤外線）と赤色のバンドを比較することで計算されます。
植生はより多くのNIR（近赤外線）を反射し、より多くの可視性がある赤色光の光を吸収します。なので、NDVI（正規化植生指標）は、植生の密度と健康状態に対する深い理解を提供します。

The formula to calculate NDVI is:

NDVIの式は以下の通り：

• NIR (Near Infrared) corresponds to Band 8 for Sentinel-2 imagery.
• Red corresponds to Band 4 for Sentinel-2 imagery.

• NIR（近赤外線） Sentinel-2衛星画像のバンド8に対応する
• 赤色 Sentinel-2衛星画像のバンド2に対応する

NDVI values range from -1 to +1, where higher values indicate healthier and denser vegetation.

NDVIの値の幅は、-1から1で、高い値はより健康かつより蜜度の高い植生です。

Sentinel-2 Band Reference

For Sentinel-2 imagery, the bands are predefined and consistent across products. Here is a quick reference for the Sentinel-2 bands:

Sentinel-2衛星画像について、バンドはあらかじめ定義されていて一貫性があります。Sentinel-2衛星画像のバンドの早見表：

(2a) - Read the raster and get the metadata (1 point)

As the first step read the raster file from data directory data/sentinel2_Nuuksio_small.tif as sentinel_raster. Provide the following information from the raster (use the same variable names and format):

ラスタファイルを、データフォルダdata/sentinel2_Nuuksio_small.tifから、sentinel_rasterとして読み込む最初のステップとして、読み込んだラスタファイルから以下の情報を取得します。（同じ変数名とフォーマットを使う。）

data_crs = "the CRS of our raster data"
resolution = "as string in this format X*Y"
total_pixels = "total number of pixels across all bands as integer"

Remember to print the values.
値を表示することを忘れずに。

import rioxarray

sentinel_raster_file = 'data/sentinel2_Nuuksio_small.tif'

sentinel_raster = rioxarray.open_rasterio(sentinel_raster_file)

data_crs = sentinel_raster.rio.crs
resolution = sentinel_raster.rio.resolution()[0] * sentinel_raster.rio.resolution()[1]
total_pixels =  len(sentinel_raster.band)  * len(sentinel_raster.x) * len(sentinel_raster.y) 

print(f"data_crs={data_crs}")
print(f"resolution={resolution}")
print(f"total_pixels={total_pixels}")

(2b) - Map NIR and RED bands (2 points)

Visualize the Near Infrared (NIR) and Red bands as individual raster plots (use subplots). You can decide what color you want to us efor your visulizations.

NIR（近赤外線）と赤色を個々のラスタのプロットとして可視化します。（サブプロットを使います。）
可視化のための使う色は自分で選ぶことができます。

import matplotlib.pyplot as plt
band_red = 4
band_NIR = 8

_, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 10))

sentinel_raster.sel(band=band_red).plot(ax=ax1, cmap="terrain")
ax1.set_title("Red Band")
ax1.set_ylabel("y cord")
ax1.set_xlabel("x cord")
sentinel_raster.sel(band=band_NIR).plot(ax=ax2, cmap="terrain")
ax2.set_title("NIR Band")
ax2.set_ylabel("y cord")
ax2.set_xlabel("x cord")

(2c) - Calculate the NDVI (5 points)

Calculate the NDVI using the formula provided above. Store the resulting raster as NDVI and export it as a .tif file. Finally, plot the NDVI raster.

前述の式を使って、NDVI（正規化植生指標）を計算します。結果のラスタをNDVIとして保存し、ひとつの .tifファイルに保存します。最後に、NDVI（正規化植生指標）のラスタをプロットします。

band_red = 4
band_NIR = 8

band_red_raster = sentinel_raster.sel(band=band_red)
band_NIR_raster = sentinel_raster.sel(band=band_NIR)
ndvi_raster = (band_NIR_raster - band_red_raster) / (band_NIR_raster + band_red_raster)

#  (nir_data + red_data) が 0 になる要素だけNaN に置き換える
ndvi_raster = ndvi_raster.where((band_NIR_raster + band_red_raster) != 0)

ndvi_raster.plot(cmap="terrain")

(2d) - Reclassify the NDVI raster (4 points)

In this task, you will reclassify the NDVI values into distinct vegetation categories and visualize the result.

このタスクでは、NDVI（正規化植生指標）の値を、独立した植生のカテゴリーに再分類し、結果を可視化します。

1. Reclassify the NDVI values:
   • Create the following NDVI classification:
     • NDVI < 0.2: Non-Vegetated
     • 0.2 ≤ NDVI < 0.5: Sparse Vegetation
     • 0.5 ≤ NDVI: Dense Vegetation
   • Use appropriate values to represent each class (1 for Non-Vegetated, 2 for Sparse Vegetation, 3 for Dense Vegetation).
2. Visualize the Reclassified NDVI:
   • Plot the reclassified NDVI using a discrete color map where each class has a distinct color.
   • Use a suitable color scheme to clearly differentiate between the classes (e.g., brown for Non-Vegetated, yellow for Sparse Vegetation, green for Dense Vegetation).

1.NDVI（正規化植生指標）の値を再分類：

• 以下のNDVI（正規化植生指標）再分類を実施：
  • NDVI < 0.2: 非植生
  • 0.2 ≤ NDVI < 0.5: まばらな植生
  • 0.5 ≤ NDVI: 密度の高い植生
• 各々の分類（1は非植生、２はまばらな植生、３は密度の高い植生）を表現するために適切な値を使うこと

これについては、xarrayライブラリのwhere関数を用います。

2.再分類したNDVI（正規化植生指標）を可視化する：

• 各々の分類が違う色を持つカラーマップを使って、再分類されたNDVI（正規化植生指標）をプロットする
• 各分類を明確に差別化するために、適切な色スキーマを使う（例えば、茶色は非植生、黄色はまばらな植生、緑は密度の高い植生）

こちらについては、Lesson 7 Processing and Analysis of Raster Dataにある方法を用います。

import xarray as xr
from matplotlib.colors import ListedColormap

NON_VEGETATED = 0
SPARSE_VEGETATION = 1
DENCE_VEGETATION = 2

classified_ndvi_raste = xr.where(
              ndvi_raster < 0.2,
              NON_VEGETATED,
              xr.where(
                  ndvi_raster < 0.5,
                  SPARSE_VEGETATION,
                  DENCE_VEGETATION
              )
)

cmap = ListedColormap(["brown", "yellow", "green"])

plt = classified_ndvi_raste.plot(cmap=cmap)
plt.colorbar.set_ticks([NON_VEGETATED, SPARSE_VEGETATION, DENCE_VEGETATION])
plt.colorbar.set_ticklabels(["Non-Vegetated","Sparse Vegetation","Dense Vegetation"])

次に、ndvi.tifファイルとして、dataディレクトリに保存します。

classified_ndvi_raste.rio.to_raster("data/ndvi.tif")

確認として、保存したファイルを読み込んでみましょう。

rioxarray.open_rasterio('data/ndvi.tif').plot()

良さそうですね。

Answer these questions:

• Q1: How does the distribution of vegetation classes compare across the area?
• Q2: Are there more areas with sparse or dense vegetation?

• Q1： エリア全体で、各植生クラスの分布は、どう比較されているか？
• Q2： まばらな植生の範囲、密度の高い植生、どちらのエリアが多いか？

-- Your answers
あなたの答え
Q1: 密度の高い植生のエリアが多い、非植生のエリアは前提的
Q2: 密度の高い植生のエリアのほうが多い