Cartographie des zones brûlées en temps quasi réel à l’aide du datacube Digital Earth Africa

• Produits utilisés : s2_l2a, ls8_sr, ls9_sr

Keywords: data used; sentinel-2, data used; landsat, band index; NBR, fire mapping

Aperçu

Ce carnet montre un exemple d’utilisation des images en temps quasi réel disponibles gratuitement sur Digital Earth Africa pour les utilisateurs qui suivent les récents incendies. Ce carnet décrit un flux de travail pour sélectionner les images « en temps quasi réel » et de référence les plus adaptées à la comparaison à partir des capteurs Landsat 8 et 9 et Sentinel-2A et Sentinel-2B. Un seuillage est appliqué pour détecter les zones brûlées non vérifiées, avant de créer un fichier vectoriel polygonal exportable.

Que signifie la cartographie des zones brûlées « en temps quasi réel » ?

Dans ce contexte, les termes « en temps quasi réel (NRT) » ou « cartographie rapide » font référence aux données ou images satellitaires des satellites Sentinel-2 et Landsat qui seront traitées dès que possible après réception des données par Digital Earth Africa. Il peut s’écouler quelques jours entre la capture d’une image par un capteur et sa mise à disposition pour analyse dans l’environnement Sandbox de Digital Earth Africa. Bien que Digital Earth Africa fasse tout son possible pour rendre les images disponibles dès que possible, des difficultés imprévues peuvent entraîner des retards dans la disponibilité des images.

Quand utiliser ce carnet

Ce carnet a pour but de cartographier l’étendue des incendies en Afrique au cours des quinze derniers jours, en fonction de la disponibilité d’images appropriées. Les résultats de ce carnet ne mesurent pas la gravité des incendies.

Pour les utilisateurs intéressés par la cartographie des incendies historiques, veuillez plutôt utiliser le carnet de données « Cartographie des zones brûlées à l’aide du carnet de données Sentinel-2 <Burnt_area_mapping.ipynb> ».

Taux de combustion normalisé (NBR) et taux de combustion normalisé delta (dNBR)

Le taux de brûlage normalisé (NBR) est un indice de gravité des incendies (FSI) qui utilise les différences dans la façon dont la végétation verte saine et la végétation brûlée réfléchissent la lumière pour détecter les pixels brûlés dans l’imagerie multispectrale. L’indice NBR nécessite des signaux des parties NIR (proche infrarouge) et SWIR (infrarouge à ondes courtes) du spectre électromagnétique et est défini ci-dessous.

\begin{equation} \text{NBR} = \frac{(\text{NIR} - \text{SWIR})}{(\text{NIR} + \text{SWIR})} \end{equation}

La comparaison des valeurs NBR les plus récentes avec une valeur de référence ou une valeur NBR antérieure (c’est-à-dire dNBR) peut aider à éliminer le bruit et à isoler les changements environnementaux dans un flux de travail EO. Le changement de NBR est appelé delta NBR (dNBR) et est défini comme suit :

\begin{equation} \text{dNBR} = \text{NBR}_{\text{ligne de base}} - \text{NBR}_{\text{après incendie}} \end{equation}

Vous trouverez plus d’informations sur le NBR et le dNBR dans le carnet de données « Cartographie des zones brûlées à l’aide de Sentinel-2 <Burnt_area_mapping.ipynb> » à la place.

Taux de combustion relativisé (RBR)

Le rapport de brûlage relativisé (RBR) est une variante du rapport de brûlage normalisé delta relativisé (RdNBR) développé par Parks et al., 2014 <https://www.mdpi.com/2072-4292/6/3/1827> qui résout certains des problèmes numériques de l’algorithme RdNBR original. Comme le RdNBR, le RBR vise à améliorer la cartographie des zones brûlées sur les zones brûlées qui présentaient de faibles niveaux de végétation avant l’incendie en prenant en compte la mesure de référence du NBR pour éviter que ces zones ne soient exclues.

\begin{equation} \text{RBR} = \frac{dNBR}{(\text{NBR}_{\text{ligne de base}} + 1,001)}\end{equation}

Plus tard dans cette analyse, vous souhaiterez peut-être utiliser l’indice RBR au lieu de l’indice dNBR si votre zone d’intérêt (AOI) comporte des quantités considérables de zones brûlées qui présentaient à l’origine de faibles niveaux de végétation. Par exemple, l’indice RBR peut être plus efficace pour cartographier les incendies sur des AOI comportant à la fois des prairies stériles et des canopées denses que la méthode dNBR.

Description

Ce cahier contient les étapes suivantes :

1. Démarrer et définir une zone d’intérêt (AOI)

2. Définir des plages de dates appropriées pour les images de base et en temps quasi réel

3. Charger, visualiser et sélectionner une image en temps quasi réel

4. Charger et sélectionner des images de base

5. Calculer les tableaux de données NBR et dNBR et effectuer un post-traitement facultatif

6. Convertir des données raster en vecteur et exporter des produits

Commencer

Pour exécuter cette analyse, exécutez toutes les cellules du bloc-notes, en commençant par la cellule « Charger les packages ».

Importez les packages Python utilisés pour l’analyse.

import math
import sys

import datacube
import numpy as np
import pandas as pd
import xarray as xr
import geopandas as gpd
import matplotlib.pyplot as plt

from scipy import ndimage
from datacube.utils import cog
from skimage import morphology
from datetime import datetime, timedelta
from odc.geo.geom import Geometry

from deafrica_tools.datahandling import load_ard
from deafrica_tools.plotting import display_map, rgb
from deafrica_tools.areaofinterest import define_area
from deafrica_tools.bandindices import calculate_indices
from deafrica_tools.spatial import xr_rasterize, xr_vectorize

Connectez-vous au datacube pour que nous puissions accéder aux données DEA. Le paramètre « app » est un nom unique pour l’analyse qui est basé sur le nom du fichier du notebook.

dc = datacube.Datacube(app="Burnt_area_mapping_near_realtime")

La cellule suivante définit les paramètres qui définissent la zone d’intérêt sur laquelle effectuer l’analyse. #### Sélectionner l’emplacement Pour définir la zone d’intérêt, deux méthodes sont disponibles :

1. En spécifiant la latitude, la longitude et la zone tampon, ou en séparant les zones tampons de latitude et de longitude, cette méthode vous permet de définir une zone d’intérêt autour d’un point central. Vous pouvez saisir la latitude centrale, la longitude centrale et une valeur de zone tampon en degrés pour créer une zone carrée autour du point central. Par exemple, « lat = 10.338 », « lon = -1.055 » et « buffer = 0.1 » sélectionneront une zone avec un rayon de 0,1 degré carré autour du point avec les coordonnées « (10.338, -1.055) ».

   Vous pouvez également fournir des valeurs de tampon distinctes pour la latitude et la longitude d’une zone rectangulaire. Par exemple, « lat = 10.338 », « lon = -1.055 » et « lat_buffer = 0.1 » et « lon_buffer = 0.08 » sélectionneront une zone rectangulaire s’étendant sur 0,1 degré au nord et au sud et 0,08 degré à l’est et à l’ouest à partir du point « (10.338, -1.055) ».

   Pour des temps de chargement raisonnables, définissez la mémoire tampon sur « 0,1 » ou moins.

2. By uploading a polygon as a GeoJSON or Esri Shapefile. If you choose this option, you will need to upload the geojson or ESRI shapefile into the Sandbox using Upload Files button in the top left corner of the Jupyter Notebook interface. ESRI shapefiles must be uploaded with all the related files (.cpg, .dbf, .shp, .shx). Once uploaded, you can use the shapefile or geojson to define the area of interest. Remember to update the code to call the file you have uploaded.

Pour utiliser l’une de ces méthodes, vous pouvez décommenter la ligne de code concernée et commenter l’autre. Pour commenter une ligne, ajoutez le symbole "#" avant le code que vous souhaitez commenter. Par défaut, la première option qui définit l’emplacement à l’aide de la latitude, de la longitude et du tampon est utilisée.

Si vous exécutez le bloc-notes pour la première fois, conservez les paramètres par défaut ci-dessous. Cela montrera comment fonctionne l’analyse et fournira des résultats significatifs. L’exemple porte sur un incendie dans une zone de la région forestière du Cap Sud, près de George, en Afrique du Sud. Fin octobre et début novembre 2018, des conditions météorologiques relativement sèches ont prévalu, créant des conditions favorables aux incendies de forêt qui ont détruit de vastes zones de forêts de plantation. Deux incendies de forêt se sont produits successivement : le premier a éclaté en raison d’un incendie d’origine humaine le 24 octobre, et le second a été déclenché par la foudre le 29 octobre.

Pour exécuter le notebook pour une zone différente, assurez-vous que les données Sentinel-2 ou Landsat 8 ou 9 sont disponibles pour la zone choisie à l’aide de DEAfrica Explorer.

# Method 1: Specify the latitude, longitude, and buffer)
aoi = define_area(lat=-33.8875, lon=22.569403, buffer=0.1)

# Method 2: Use a polygon as a GeoJSON or Esri Shapefile.
# aoi = define_area(vector_path='aoi.shp')

#Create a geopolygon and geodataframe of the area of interest
geopolygon = Geometry(aoi["features"][0]["geometry"], crs="epsg:4326")
geopolygon_gdf = gpd.GeoDataFrame(geometry=[geopolygon], crs=geopolygon.crs)

# Get the latitude and longitude range of the geopolygon
lat_range = (geopolygon_gdf.total_bounds[1], geopolygon_gdf.total_bounds[3])
lon_range = (geopolygon_gdf.total_bounds[0], geopolygon_gdf.total_bounds[2])

Afficher l’emplacement sélectionné

La cellule suivante affichera la zone sélectionnée sur une carte interactive. N’hésitez pas à zoomer et dézoomer pour mieux comprendre la zone que vous allez analyser. En cliquant sur n’importe quel point de la carte, vous découvrirez les coordonnées de latitude et de longitude de ce point.

display_map(x=lon_range, y=lat_range)

Make this Notebook Trusted to load map: File -> Trust Notebook

Par défaut, les images des 14 jours précédents à partir de la date actuelle sont extraites du datacube et visualisées pour une sélection manuelle. Si une plage de dates NRT différente est souhaitée, saisissez-la pour la variable « nrt_time_delta » au format « DD ».

Pour la sélection d’une image de référence, la valeur par défaut consiste à extraire et à visualiser les images acquises au cours de la période comprise entre 14 jours et deux mois avant la date actuelle. Si une plage de dates différente est préférée, ajustez le « delta_de_temps_de_baseline » ci-dessous au format « DD ».

Les plages de dates calculées ci-dessous constituent un bon point de départ, mais l’AOI que vous avez choisi peut nécessiter des plages de dates plus grandes ou plus petites.

Afin de démontrer la capacité de ce bloc-notes lors d’un incendie, une date statique a été saisie pour la variable « nrt_date_end ». Pour exécuter ce bloc-notes pour un incendie en cours, exécutez plutôt la deuxième ligne de code notée ci-dessous pour définir « nrt_date_end » sur la date du jour.

# Near Real-time event date '2018-11-12' is used to demonstrate the notebooks capabilities.
# For use mapping a current event, please move the '#' from the second to the first line below.
nrt_date_end = "2018-11-12"
# nrt_date_end = datetime.today().strftime("%Y-%m-%d")

# Define the date ranges for NRT and baseline imagery requests. The default is the preceding 14 days for
# NRT imagery, then the two months prior to that for baseline imagery.

nrt_time_delta = 14
baseline_time_delta = 42

# Calculate the beginning of the 'nrt' date range.
nrt_date_start = datetime.strftime(
    ((datetime.strptime(nrt_date_end, "%Y-%m-%d")) - timedelta(days=nrt_time_delta)),
    "%Y-%m-%d",
)
# Date range from which to extract and visualise potential baseline images. Can be manually changed by enterring a date in 'YYYY-MM-DD' format.
baseline_start = datetime.strftime(
    (
        (datetime.strptime(nrt_date_start, "%Y-%m-%d"))
        - timedelta(days=baseline_time_delta)
    ),
    "%Y-%m-%d",
)
baseline_end = datetime.strftime(
    ((datetime.strptime(nrt_date_end, "%Y-%m-%d")) - timedelta(days=nrt_time_delta)),
    "%Y-%m-%d",
)
# Print Date ranges
print(
    f"Potential Near Real-time Images will be extracted between {nrt_date_start} and {nrt_date_end}"
)
print(
    f"Potential baseline Images will be extracted between {baseline_start} and {baseline_end}"
)

Potential Near Real-time Images will be extracted between 2018-10-29 and 2018-11-12
Potential baseline Images will be extracted between 2018-09-17 and 2018-10-29

La disponibilité d’images récentes et sans nuages est le facteur limitant le plus important pour la réalisation d’une cartographie des zones brûlées en temps quasi réel lorsqu’un incendie se produit. Par conséquent, nous examinerons les images récentes des collections d’images Sentinel-2 et Landsat pour maximiser les chances de trouver une image appropriée.

Pour réduire l’impact de l’influence des capteurs croisés sur notre analyse, nous chargerons uniquement les images de base du capteur à partir duquel notre image en temps quasi réel est sélectionnée.

Définir les paramètres de charge

NB : les variables « min_gooddata_nrt » et « min_gooddata_base » sont utilisées pour filtrer les images de mauvaise qualité pour l’AOI choisi dans les extraits d’images en temps quasi réel et de base respectivement. Les valeurs attribuées à ces variables définissent le seuil du nombre de pixels « bons » (c’est-à-dire sans nuages) requis pour que l’image soit extraite.

Une valeur « min_gooddata » de 0,9 garantit que seules les images contenant plus de 90 % de pixels sans nuages sur l’AOI sont extraites. Pour augmenter le nombre d’images extraites, ces valeurs peuvent être réduites pour permettre l’extraction d’images contenant plus de pixels nuageux.

# Set spatial, spectral, and quality parameters for the imgery extracrt from the datacube.

# Setting the min_gooddata value lower for the near real-time extract increases the number of images to chose from
min_gooddata_nrt = 0.6
min_gooddata_base = 0.8

# DIfferent measurements are required for each sensor due to changes in band nomenclature
measurements_s2 = [
    "blue",
    "green",
    "red",
    "nir_1",
    "swir_2",
]
measurements_ls = [
    "blue",
    "green",
    "red",
    "nir",
    "swir_2",
]

# Define the resolution tobe used for each sensor. The Sentinel-2 resolution can be reduced to 10m if desired.
s2_res = (-20, 20)
ls_res = (-30, 30)

# Define the coordinate system
output_crs = "EPSG:6933"

# Create a query object for the universal parameters
query = {
    "x": lon_range,
    "y": lat_range,
    "output_crs": output_crs,
    "group_by": "solar_day",
}

Charger des images de la collection Sentinel-2

# Load imagery from the Sentinel-2 Collection. Note the s2cloudless mask is applied instead of the standard fmask.
nrt_s2 = load_ard(
    dc=dc,
    products="s2_l2a",
    time=(nrt_date_start, nrt_date_end),
    measurements=measurements_s2,
    min_gooddata=min_gooddata_nrt,
    resolution=s2_res,
    **query
)

#Assign a custom attribute 'sensor' to the Sentinel-2 dataset for plotting purposes.
nrt_s2 = nrt_s2.assign_attrs({"sensor":"s2"})

Using pixel quality parameters for Sentinel 2
Finding datasets
    s2_l2a
Counting good quality pixels for each time step

/opt/venv/lib/python3.12/site-packages/rasterio/warp.py:387: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  dest = _reproject(

Filtering to 5 out of 6 time steps with at least 60.0% good quality pixels
Applying pixel quality/cloud mask
Loading 5 time steps

Charger des images de la collection Landsat

# Load imagery from the Landsat Collection.
nrt_ls = load_ard(
    dc=dc,
    products=["ls8_sr", "ls9_sr"],
    time=(nrt_date_start, nrt_date_end),
    measurements=measurements_ls,
    min_gooddata=min_gooddata_nrt,
    resolution=ls_res,
    **query
)

#Assign a custom attribute 'sensor' to the Landsat dataset for plotting purposes.
nrt_ls = nrt_ls.assign_attrs({"sensor":"ls"})

Using pixel quality parameters for USGS Collection 2
Finding datasets
    ls8_sr
    ls9_sr
Counting good quality pixels for each time step
Filtering to 1 out of 2 time steps with at least 60.0% good quality pixels
Applying pixel quality/cloud mask
Re-scaling Landsat C2 data
Loading 1 time steps

Visualisation d’images en temps quasi réel

Les images des collections Sentinel-2 et Landsat qui répondent aux exigences spatiales, temporelles et de qualité que nous avons définies précédemment seront visualisées ci-dessous. Notez la date, le numéro d’index et le capteur de chaque image située dans le titre de la sous-parcelle.

Dans l’étape suivante, nous devrons définir les variables « nrt_img_index » et « nrt_sensor » sur l’index des images et le capteur choisis respectivement.

Il n’est pas rare qu’aucune image Landsat appropriée n’ait été acquise au cours des 14 jours précédents et qu’elle réponde aux exigences en temps quasi réel. Si cela se produit, le notebook continuera à fonctionner sans problème.

Définir la fonction « plot_variable_images »

La fonction plot_variable_images trace dynamiquement des images à partir d’un ensemble de données xarray contenant des bandes RVB, avec des titres affichant la date et l’index, et éventuellement le type de capteur s’il est spécifié. Elle vérifie d’abord que l’entrée est un ensemble de données xarray valide et calcule le nombre d’images. Si aucune image n’est présente, elle génère une erreur. La fonction détermine le nombre nécessaire de lignes de sous-tracés, crée les sous-tracés et parcourt chaque image pour la tracer sous forme d’image RVB, garantissant ainsi une mise à l’échelle des couleurs robuste.

def plot_variable_images(img_collection):
    # Check that img_collection is a xarray dataset
    if not isinstance(img_collection, xr.Dataset):
        raise TypeError("img_collection must be a xarray dataset.")

    # Calculate number of images in `img_collection`
    plot_count = img_collection.dims["time"]

    # Check if dataset has 0 images
    if plot_count == 0:
        if hasattr(img_collection, "sensor"):
            raise ValueError("The {} dataset has no images to display for the "
            "given query parameters".format(img_collection.sensor))
        else:
            raise ValueError("The supplied xarray dataset has no images to "
            "display for the given query parameters")

    # Divide the number of images by 2 rounding up to calculate the
    # number of rows for the below figure are needed
    plot_rows = math.ceil(plot_count / 2)

    # Construct a figure to visualise the imagery
    f, axarr = plt.subplots(plot_rows, 2, figsize=(10, plot_rows * 4.5),
                            squeeze=False)

    # Flatten the subplots so they can easily be enumerated through
    axarr = axarr.flatten()

    # Iterate through each image in the dataset and plot
    # each image as a RGB on a subplot
    for t in range(plot_count):
        rgb(
            img_collection.isel(time=t),
            bands=["red", "green", "blue"],
            ax=axarr[t],
            robust=True,
        )
        # Test to see if the dataset has a custom 'sensor' attribute.
        # If so, include the string in each subplot title.
        if hasattr(img_collection, "sensor"):
            title = (
                str(img_collection.time[t].values)[:10]
                + "  ||  Index: "
                + str(t)
                + "  ||  Sensor: "
                + img_collection.sensor
            )
        else:
            title = (
                str(img_collection.time[t].values)[:10]
                + "  ||  Index: "
                + str(t)
            )
        # Set subplot title, axis label, and shrink offset text
        axarr[t].set_title(title)
        axarr[t].set_xlabel("X coordinate")
        axarr[t].set_ylabel("Y coordinate")
        axarr[t].yaxis.offsetText.set_fontsize(6)
        axarr[t].xaxis.offsetText.set_fontsize(6)

    # Adjust padding arround subplots to prevent overlapping elements
    plt.tight_layout()

    # Remove the last subplot if an odd number of images are being displayed
    if plot_count % 2 != 0:
        f.delaxes(axarr[plot_count])

#Plot imagery from the Sentinel-2 satellites by running the `plot_variable_images` tool on the `nrt_s2` xarray dataset.
plot_variable_images(nrt_s2)

#Plot suitable imagery from the Landsat satellites by running the `plot_variable_images` tool on the `nrt_ls` xarray dataset.
if nrt_ls.time.size == 0:
    print('No suitable Landsat imagery')
else:
    plot_variable_images(nrt_ls)

En examinant les images en temps quasi réel (NRT) disponibles à la fois sur Sentinel-2 et sur Landsat, l’image Landsat semble être l’option la plus adaptée pour une analyse plus approfondie. Bien qu’il existe une seule image Sentinel-2 (Index : 3) avec une résolution spatiale plus élevée, elle présente plus de lacunes que l’image Landsat disponible.

Pour le processus de sélection de l’image de référence qui se déroule quelques cellules plus loin dans le code, les images du capteur de l’image NRT sélectionnée seront chargées. Si une image Sentinel-2 NRT est sélectionnée, les images Sentinel-2 seront chargées pour la sélection de l’image de référence. De même, si une image Landsat NRT est choisie, les images Landsat seront chargées pour la sélection. Cela garantit la continuité et la cohérence de l’ensemble de données, facilitant une comparaison et une analyse précises entre les périodes NRT et de référence.

Sélection d’images en temps quasi réel

À partir des images ci-dessus, notez la valeur d’index de la scène la plus appropriée et définissez la variable « nrt_img_index » ci-dessous sur cette valeur. Par exemple, si la première image affichée est la plus désirable, définissez « nrt_img_index = 0 ».

De plus, définissez la variable rt_sensor sur le capteur d’où provient l’image en temps quasi réel choisie, s2 pour Sentinel-2 ou ls pour Landsat. Seules les images de la collection sélectionnée seront affichées pour la visualisation de l’image de base ci-dessous afin de réduire l’influence des différents capteurs sur l’analyse.

# Index selection for the Near Real-time image chosen above. The default value of -1 selects the most recent image regardless of suitability.
nrt_img_index = 0
# Set the below variable to either 's2' or 'ls' depending on the sensor your chosen image is from.
nrt_sensor = "ls"

# Assign the selected image to the rt_img variable, and identify which sensor and measurements to use to find a baseline image.
if nrt_sensor == "s2":
    nrt_img = nrt_s2.isel(time=nrt_img_index)
    baseline_products = "s2_l2a"
    baseline_measurements = measurements_s2
    index_collection = "s2"
    baseline_res = s2_res
elif nrt_sensor == "ls":
    nrt_img = nrt_ls.isel(time=nrt_img_index)
    baseline_products = ["ls8_sr", "ls9_sr"]
    baseline_measurements = measurements_ls
    index_collection = "c2"
    baseline_res = ls_res
else:
    print(
        "Please make sure the 'nrt_sensor' variable is correctly set to either s2 or ls"
    )

Charger et sélectionner des images de base

# Load imagery from whichever collection the NRT image was selected from.
baseline = load_ard(
    dc=dc,
    products=baseline_products,
    time=(baseline_start, baseline_end),
    measurements=baseline_measurements,
    min_gooddata=min_gooddata_base,
    resolution=baseline_res,
    **query
)

#Assign a custom attribute 'sensor' to the baseline dataset for plotting purposes.
baseline = baseline.assign_attrs({"sensor": nrt_sensor})

Using pixel quality parameters for USGS Collection 2
Finding datasets
    ls8_sr
    ls9_sr
Counting good quality pixels for each time step
Filtering to 3 out of 5 time steps with at least 80.0% good quality pixels
Applying pixel quality/cloud mask
Re-scaling Landsat C2 data
Loading 3 time steps

Visualisez les images de base extraites

#Plot imagery from the baseline dataset by running the `plot_variable_images` tool on the `baseline` xarray dataset.
plot_variable_images(baseline)

Sélection de l’image de base

Il faut choisir une image de référence parmi celles sélectionnées ci-dessus. Il est important de se rappeler que l’indice dNBR mesure les changements environnementaux associés aux changements de paysage induits par les incendies. Par conséquent, une scène similaire à l’image en temps quasi réel, exempte d’effets d’incendie, fournira le meilleur point de comparaison pour détecter les cicatrices de brûlures. Par exemple, si le paysage de l’image en temps quasi réel est sec, il est préférable de choisir une image de référence où le paysage présente un niveau de sécheresse similaire à celui d’une image plus verte.

Notez le numéro d’index de l’image de base choisie ci-dessus, puis définissez-le sur la variable baseline_img_index.

# Set the baseline_img_index variable to the index of the chosen baseline image
baseline_img_index = 2

# Extract the chosen baseline image to the baseline_img variable
baseline_img = baseline.isel(time=baseline_img_index)

Calculer les tableaux de données NBR et dNBR et effectuer un post-traitement facultatif

# Calculate NBR for the near real-time image and assign it to the rt_NBR variable
nrt_image_NBR = calculate_indices(
    nrt_img, index="NBR", collection=index_collection, drop=False
)
nrt_NBR = nrt_image_NBR.NBR

# Calculate NBR for the baseline image and assign it to the baseline_NBR variable
baseline_image_NBR = calculate_indices(
    baseline_img, index="NBR", collection=index_collection, drop=False
)
baseline_NBR = baseline_image_NBR.NBR

# Calculate dNBR (delta NBR) by differeing the two indices
dNBR = baseline_NBR - nrt_NBR

# Create a figure to visualise the selected images in true colour, as well as the NBR index.
f, axarr = plt.subplots(2, 2, figsize=(12, 10), squeeze=False, layout="constrained")

# Visualise the selected near real-time and baseline images in true colour
rgb(
    nrt_img,
    bands=["red", "green", "blue"],
    ax=axarr[0, 0],
    robust=True,
)

rgb(
    baseline_img,
    bands=["red", "green", "blue"],
    ax=axarr[1, 0],
    robust=True,
)

# Visualise the NBR index for each image
nrt_NBR.plot(cmap="RdBu", vmin=-1, vmax=1, ax=axarr[0, 1])

baseline_NBR.plot(cmap="RdBu", vmin=-1, vmax=1, ax=axarr[1, 1])

# Set subplot Titles
axarr[0, 0].set_title("Near Real-time RGB " + str(nrt_img.time.values)[:10])
axarr[0, 1].set_title("Near Real-time NBR " + str(nrt_img.time.values)[:10])
axarr[1, 0].set_title("Baseline RGB " + str(baseline_img.time.values)[:10])
axarr[1, 1].set_title("Baseline NBR " + str(baseline_img.time.values)[:10])

Text(0.5, 1.0, 'Baseline NBR 2018-10-24')

dNBR vs RBR pour la cartographie des zones brûlées

Comme nous l’avons vu précédemment, deux indices de gravité des incendies différenciés sont disponibles dans ce carnet. Essayez d’utiliser l’indice dNBR dans un premier temps. Cependant, si l’incendie dans votre zone d’intérêt se produit sur des paysages à végétation moins dense, comme des prairies, vous souhaiterez peut-être tester l’indice RBR. Dans cet exemple, nous avons utilisé l’indice RBR car l’incendie s’est produit sur des surfaces importantes de zones dominées par l’herbe.

En fonction de votre choix d’indice de gravité des incendies, définissez la variable « bam_method » sur « dNBR » ou « RBR ».

# set the bam_method variable to either 'dNBR' or 'RBR'.
bam_method = "RBR"

# Calculate the bam variable based on the chosen dFSI
if bam_method == "RBR":
    RBR = dNBR / (baseline_NBR + 1.001)
    bam = RBR

elif bam_method == "dNBR":
    bam = dNBR
else:
    print(
        "Please make sure the 'bam_method' variable is correctly set to either 'RBR' or 'dNBR'"
    )

#Rename the variable in the bam dataarray with the new dFSI
bam.name=(bam_method)

# Create a figure to plot the chosen fire severity index
f, axarr = plt.subplots(
    1, 2, figsize=(15, 10), squeeze=False, gridspec_kw={"width_ratios": [1, 5]}
)

# Calculate and round the dNBR dataarray value range to set determine the plots colourmap range
bam_NBR_min = round(float(bam.quantile(0.005)), 1)
bam_NBR_max = round(float(bam.quantile(0.995)), 1)

# PLot the dNBR dataarray on the second subplot of the above figure
bam.plot(cmap="RdBu_r", vmin=bam_NBR_min, vmax=bam_NBR_max, ax=axarr[(0, 1)])

# Plot a histogram of dNBR values in the first figure subplot.
# Calculate a colourmap from the dataarray plot by iterating through individual histogram patches
cm = plt.colormaps.get_cmap("RdBu_r")

n, bins, patches = xr.plot.hist(
    darray=bam,
    bins=np.arange(bam_NBR_min, bam_NBR_max + 0.05, 0.05),
    align="mid",
    orientation="horizontal",
    ec="black",
    yticks=(np.arange(bam_NBR_min - 0.05, bam_NBR_max + 0.05, step=0.05)),
    ax=axarr[(0, 0)],
)

# Match the colour scale of the histogram to that used in the map plot.
bin_centers = 0.5 * (bins[:-1] + bins[1:])
col = bin_centers - min(bin_centers)
col /= max(col)
for c, p in zip(col, patches):
    plt.setp(p, "facecolor", cm(c))

# Set titles for each subplot
axarr[0, 0].set_title(bam_method + " Histogram")
axarr[0, 1].set_title(
    bam_method
    + " measured between "
    + str(baseline_img.time.values)[:10]
    + " - "
    + str(nrt_img.time.values)[:10]
)

# Set the x-axis label and number of x-axis ticks for the histogram plot
axarr[0, 0].set_xlabel(bam_method + " count")
axarr[0, 0].xaxis.set_major_locator(plt.MaxNLocator(3))

Définition de la valeur seuil pour identifier les zones brûlées

Il faut choisir une valeur pour délimiter les zones brûlées et non brûlées. Cette valeur peut varier en fonction de la structure de la végétation et des conditions environnementales de la zone cartographiée. Une valeur de 0,3 peut être utilisée comme point de départ. Cependant, il est utile d’interpréter l’histogramme ci-dessus et le tracé dNBR pour vous aider à déterminer la valeur la plus adaptée à votre zone d’intérêt. La sélection d’une valeur seuil implique un compromis entre une valeur si basse que les faux retours sont inclus et une valeur suffisamment élevée pour exclure les zones brûlées de faible gravité et moins végétalisées.

Définissez la valeur de seuil choisie sur la variable « seuil » ci-dessous.

# Set threshold value. Fire serverity index values below this value will be masked out.
threshold = 0.30

# Apply threshold to the dNBR dataset to create the `burnt` dataset of burnt pixels
burnt = bam > threshold

# Mask real-time true colour image with the above `burnt` mask to show what has been captured as burnt area.
masked = bam.where(burnt == 1)
bands = ["red", "green", "blue"]
rgb(nrt_img.where(burnt == 1), bands=bands)

Optionnel : post-traitement morphologique

Le résultat de notre analyse peut avoir généré une quantité inacceptable de pixels isolés et d’autres retours bruyants. Nous pouvons appliquer des opérations morphologiques au dataarray binaire xr pour supprimer cela des produits et simplifier la géométrie des polygones de sortie. Cette étape est facultative, mais peut être utile lorsque l’analyse produit des sorties bruyantes et « mouchetées ».

Les opérations de fermeture, d’érosion et de dilatation sont décrites ci-dessous. La fermeture est une opération composée d’une dilatation suivie d’une érosion pour supprimer les petits trous. L’érosion rétrécit les pixels de l’image, ce qui permet de supprimer les retours bruyants et isolés. La dilatation étend les pixels de l’image pour supprimer les petits trous et simplifier la géométrie des polygones dérivés.

Ces opérations sont réalisées à l’aide d’un élément de structuration de disque. Une taille par défaut de 3 pixels est utilisée ci-dessous, mais cette valeur peut être modifiée pour faire varier le degré de post-traitement appliqué.

# Define the size of the disk structuring element, measured in number of pixels.
# The default value is 5.
disk_size = 5

# Perform the Closing, Erosion, and Dilation operations to the burnt dataarray
dilated_data = xr.DataArray(
    morphology.binary_closing(burnt, morphology.disk(disk_size)), coords=burnt.coords
)
erroded_data = xr.DataArray(
    morphology.erosion(dilated_data, morphology.disk(disk_size)), coords=burnt.coords
)
dilated_data = xr.DataArray(
    ndimage.binary_dilation(erroded_data, morphology.disk(disk_size)),
    coords=burnt.coords,
)

# Save the results to the burnt variable
burnt = dilated_data

# Visualise the post-processed dataarray to show the final burnt area
rgb(nrt_img.where(dilated_data == 1), bands=bands)

Convertir des données raster en vecteur et exporter des produits

Exporter des rasters

Trois rasters sont exportés ci-dessous :

1. An RGB geotiff of the Near Real-time imagery

2. Un géotiff RVB de l’imagerie de base

3. Un géotiff monobande non masqué de l’indice delta NBR choisi - Ceux qui souhaitent découper ces images sur la zone brûlée délimitée peuvent utiliser le fichier de formes ci-dessous.

# Set a name for the AOI (used to name exported products)
area_name = "Example"

# Write near real-time fire image to multi-band GeoTIFF
cog.write_cog(
    geo_im=nrt_img.to_array(),
    fname=f"{area_name}_{str(nrt_img.time.values)[:10]}_near_realtime_image.tif",
    overwrite=True,
)

# Write baseline reference image to multi-band GeoTIFF
cog.write_cog(
    geo_im=baseline_img.to_array(),
    fname=f"{area_name}_{str(baseline_img.time.values)[:10]}_baseline_image.tif",
    overwrite=True,
)

# Turn delta NBR into a Xarray Dataset for export to GeoTIFF
cog.write_cog(geo_im=bam, fname=f"{area_name}_{str(nrt_img.time.values)[:10]}_{bam_method}.tif", overwrite=True)

PosixPath('Example_2018-11-09_RBR.tif')

Conversion des données raster bam au format vectoriel avant l’exportation vers shapefile

Nous allons convertir le tableau de données raster au format vectoriel à l’aide de l’outil xr_vectorize et l’exporter vers un fichier de formes. Le fichier de formes sera visualisé ci-dessous.

# Convert the burnt area from raster to vector format
gdf = xr_vectorize(
    da=burnt, mask=burnt.values == 1, output_path=f"{area_name}_{str(nrt_img.time.values)[:10]}_{bam_method}.shp"
)

# Plot the vector data overlying a transparent NRT RGB Imagery
f, ax = plt.subplots(figsize=(8, 7))

gdf.plot(edgecolor="black", color="orange", ax=ax)

rgb(nrt_img, bands=["red", "green", "blue"], ax=ax, alpha=0.5)

ax.set_title(
    f"Vectorised Unverified Burnt Area Polygon for \n{area_name} measured on "
    + str(nrt_img.time.values)[:10]
)

Exporting vector data to Example_2018-11-09_RBR.shp

Text(0.5, 1.0, 'Vectorised Unverified Burnt Area Polygon for \nExample measured on 2018-11-09')

---

Informations Complémentaires

Licence Le code de ce notebook est sous licence Apache, version 2.0 <https://www.apache.org/licenses/LICENSE-2.0>`__.

Les données de Digital Earth Africa sont sous licence Creative Commons Attribution 4.0 <https://creativecommons.org/licenses/by/4.0/>.

Contact Si vous avez besoin d’aide, veuillez poster une question sur le canal Slack DE Africa <https://digitalearthafrica.slack.com/> ou sur le GIS Stack Exchange <https://gis.stackexchange.com/questions/ask?tags=open-data-cube> en utilisant la balise open-data-cube (vous pouvez consulter les questions posées précédemment `ici <https://gis.stackexchange.com/questions/tagged/open-data-

Si vous souhaitez signaler un problème avec ce notebook, vous pouvez en déposer un sur Github.

Version de Datacube compatible

print(datacube.__version__)

1.8.20

Dernier test :

from datetime import datetime
datetime.today().strftime('%Y-%m-%d')

'2025-01-15'