Analiza i monitorowanie powodzi przy użyciu języka Python i zdjęć satelitarnych Sentinel-2 na NSIS

Skoncentrujemy się na uzyskaniu obrazów obszarów dotkniętych powodzią w środkowej Grecji we wrześniu 2023 roku.

Wymagania wstępne

Nr 1 Hosting

Jest wymagane konto hostingowe NSIS z dostępem do interfejsu Horizon: https://horizon.cloudferro.com.

Nr 2 Instalacja Jupyter Notebook

Kod przedstawiony w tym artykule działa w środowisku Jupyter Notebook i można go zainstalować na wybranej platformie, postępując zgodnie z informacjami podanymi na oficjalnej stronie instalacji Jupyter Notebook.

Jednym ze szczególnych sposobów instalacji (choć w żadnym wypadku nie jest to wymagane do śledzenia tego artykułu) jest zainstalowanie aplikacji na klastrze Kubernetes. W takim przypadku będzie przydatny artykuł:

Instalacja JupyterHub na klastrze Magnum Kubernetes w NSIS Cloud.

Krok 1 Wybór odpowiednich obrazów

Pierwszym krokiem po pobraniu produktów jest wybranie odpowiednich obrazów. Do analizy MNDWI będziemy potrzebować obrazów satelitarnych w widzialnym zakresie widma, zazwyczaj o

długości fali od 490 do 610 nanometrów (ZIELONY, PASMO 3), jak również w

krótkofalowej podczerwieni (SWIR, PASMO 11).

Ponadto użyjemy Scene Classification Layer (SCL), aby zamaskować pokrycie chmurami, które może mieć wpływ na końcowy wynik analizy. Następnie otworzymy wszystkie obrazy w notatniku w języku Python.

#import libaries
import rasterio
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap

input_band_3_b = 'band_3_before_flood.jp2'
input_band_11_b = 'band_11_before_flood.jp2'
input_scl_b = 'scl_image_before_flood.jp2'

input_band_3_f = 'band_3_during_flood.jp2'
input_band_11_f = 'band_11_during_flood.jp2'
input_scl_f = 'scl_image_during_flood.jp2'

Otwieranie obrazów rastrowych

def open(raster):
    with rasterio.open(raster) as src:
        raster = src.read(1)
    return raster

#Get CRS and transformation data
with rasterio.open(input_scl_f) as src:
    crs = src.crs
    transform = src.transform

band_green_b= open(input_band_3_b)
band_swir_b = open(input_band_11_b)
scl_b = open(input_scl_b)

band_green_f = open(input_band_3_f)
band_swir_f = open(input_band_11_f)
scl_f = open(input_scl_f)

Krok 2 Obliczanie indeksu MNDWI dla wybranych obrazów

Normalizacja wartości pikseli w obrazach rastrowych

Wartości pikseli w obrazach rastrowych muszą być znormalizowane do wartości z zakresu od 0 do 255. Następnie możemy obliczyć wartość indeksu MNDWI za pomocą wzoru:

MNDWI = (ZIELONY - SWIR) / (ZIELONY + SWIR)

def normalize(band_green, band_swir):
    #Calculate maximum values for each raster file
    max_value = band_green.max()
    max_value = band_swir.max()

    #Normalize raster data to 0-255 scale
    scaled_band_green = (band_green / max_value) * 255
    scaled_band_swir = (band_swir/ max_value) * 255

    return scaled_band_green, scaled_band_swir

#Calculate MNDWI value for given images
def mndwi(scaled_band_green, scaled_band_swir):
    mndwi = (scaled_band_green - scaled_band_swir) / (scaled_band_green + scaled_band_swir)
    return mndwi

band_green_b, band_swir_b = normalize(band_green_b, band_swir_b)
band_green_f, band_swir_f = normalize(band_green_f, band_swir_f)

mndwi_b = mndwi(band_green_b, band_swir_b)
mndwi_f = mndwi(band_green_f, band_swir_f)

Wyświetlanie wyników na obrazie matplotlib

# Define a colormap that transitions from white to light green to blue to dark navy blue
colors = [(1, 1, 1), (0.8, 1, 0.8), (0, 0, 1), (0, 0, 0.4)]  # (R, G, B)
cmap = LinearSegmentedColormap.from_list('custom_colormap', colors, N=256)

# Create a figure and axis for the MNDWI before flood (mndwi_b) raster
plt.figure(figsize=(8, 8))
plt.imshow(mndwi_b, cmap=cmap, vmin=-1, vmax=1)  # Adjust the colormap and limits as needed
plt.colorbar(label='MNDWI')
plt.title('MNDWI Before Flood')

# Create a figure and axis for the MNDWI during flood (mndwi_f) raster
plt.figure(figsize=(8, 8))
plt.imshow(mndwi_f, cmap=cmap, vmin=-1, vmax=1)  # Adjust the colormap and limits as needed
plt.colorbar(label='MNDWI')
plt.title('MNDWI During Flood')

plt.show()

To jest obraz „przed powodzią”:

a to jest „po powodzi”:

Krok 3 Identyfikacja obszarów o zwiększonej obecności wody

Następnym krokiem jest uczynienie procesu bardziej przejrzystym poprzez identyfikację obszarów o zwiększonej obecności wody. Zostanie to wykonane poprzez odjęcie wartości indeksu MNDWI przed powodzią od wartości indeksu MNDWI podczas powodzi przy jednoczesnym wykluczeniu obszarów pokrytych chmurami. Wynikiem tego procesu będzie obraz, który podkreśli regiony, w których wzrosła obecność wody.

Innym cennym produktem, który możemy uzyskać, jest obraz, który w szczególności pokazuje obszary znajdujące się pod wodą podczas powodzi. Aby to osiągnąć, wykorzystamy te same obliczenia, co w poprzednim kroku, lecz z dodatkową warstwą maskowania, aby wykluczyć obszary, które miały dodatni wskaźnik MNDWI przed powodzią. W wyniku tego zostanie utworzony obrazu pokazujący obszary, w których obecnie znajduje się woda, a w których wcześniej jej nie było.

#Create mask by SCL values
mask_f = np.logical_or(scl_f == 8, np.logical_or(scl_f == 9, scl_f == 3))
mask_b = np.logical_or(scl_b == 8, np.logical_or(scl_b == 9, scl_b == 3))
#Calculate difference between flood period and normal period
diff = mndwi_f - mndwi_b

#Calculate water tides by showing the pixcels that have positive difference values
tide = diff
tide[(diff < 0)] = 0
tide[mask_f | mask_b] = 0
tide = tide / tide.max()
#Calculate flood areas by showing the pixcels that have possitive difference values and had negative values of MNDWI before flood
flood = diff
flood[(diff < 0) | (mndwi_b > 0)] = 0
flood[mask_f | mask_b] = 0
flood = flood / flood.max()

Wyświetlanie wyników na obrazie matplotlib

# Define a colormap that goes from yellow to red with transparency for 0 values and violet for high values
colors = [(1, 1, 0, 0), (1, 0.5, 0, 1), (1, 0, 0, 1), (0.5, 0, 0.5, 1)]  # (R, G, B, Alpha)
cmap = LinearSegmentedColormap.from_list('custom_colormap', colors, N=256)

# Create a figure and axis for the Water Tides raster
plt.figure(figsize=(8, 8))
plt.imshow(tide, cmap=cmap, vmin=0, vmax=np.max(tide))  # Adjust the colormap and limits as needed
plt.colorbar(label='Tide')
plt.title('Water Tides')

# Create a figure and axis for the Flood Areas raster
plt.figure(figsize=(8, 8))
plt.imshow(flood, cmap=cmap, vmin=0, vmax=np.max(flood))  # Adjust the colormap and limits as needed
plt.colorbar(label='Flood')
plt.title('Flood Areas')

plt.show()

Obraz dla fali powodziowych:

Obraz dla zalanych terenów:

Pobierz obrazy wynikowe

Ostatnim krokiem będzie pobranie wynikowych obrazów przy użyciu biblioteki Rasterio:

def save(raster, name):
    with rasterio.open(f"{name}.tif", 'w', driver='GTiff', width=raster.shape[1], height=raster.shape[0], count=1, dtype=raster.dtype, crs=crs, transform=transform) as dst:
        dst.write(raster, 1)

save(mndwi_b, 'mndwi_before_greece')
save(mndwi_f, 'mndwi_flood_greece')
save(tide, 'tide_greece')
save(flood, 'flood_greece')