Blog SIG & Territoires

Dans deux articles précédents, nous avons exploré comment passer des images Sentinel-2 standard (10 m de résolution) à une version super-résolue à 1 m à l’aide du module S2DR3, développé pour la recherche sur les coraux de Maurice. Cette étape permet d’obtenir un niveau de détail inédit à partir de données ouvertes, en révélant des structures fines souvent invisibles à 10 m : herbiers, platiers, zones sableuses, ou encore transitions subtiles entre végétation et récif.

Traditionnellement, le calcul des indices spectraux (NDVI, NDWI, MNDWI, BSI, NDBI, etc.) se fait dans QGIS, une fois les images téléchargées et traitées. Cependant, il est souvent plus efficace de générer ces indices directement à la source, c’est-à-dire dès la production des images Sentinel-2 super-résolues.

Cette approche présente plusieurs avantages :

• Automatisation complète : les indices sont calculés dès la création du jeu de données, sans traitement manuel supplémentaire dans QGIS.
• Cohérence spatiale : les indices utilisent exactement les mêmes corrections radiométriques et la même super-résolution que les bandes d’origine.
• Gain de temps : aucun export intermédiaire ni configuration manuelle d’algorithmes raster.
• Préparation pour l’analyse : les résultats sont immédiatement exploitables dans QGIS ou tout autre outil SIG.

Le script présenté ci-dessous automatise ainsi l’ensemble du flux :

1. téléchargement et super-résolution de la scène Sentinel-2,
2. génération d’un masque de nuages intelligent (avec exclusion des zones de déferlement),
3. calcul direct des principaux indices spectraux,
4. création d’un fichier multibande prêt à l’analyse.

Cette intégration simplifie considérablement le passage de l’image brute aux cartes d’interprétation, tout en garantissant la reproductibilité du traitement scientifique.

Les indices calulés

NDVI — Normalized Difference Vegetation Index

• But : mesurer la densité et la santé de la végétation.
• Formule :

.

.

• NIR = bande proche infrarouge (B08)
• Red = bande rouge (B04)
• Valeurs : -1 à 1 ; >0.2 végétation, <0.1 sol nu ou eau.
• Usage : suivi de la végétation, estimation de la biomasse, agriculture.

NDWI — Normalized Difference Water Index

• But : détecter la présence d’eau dans le paysage.
• Formule :

.

.

• Green = bande verte (B03)
• NIR = bande proche infrarouge (B08)
• Valeurs : positives = eau, négatives = sol/végétation.
• Usage : cartographie des lacs, rivières, zones humides.

MNDWI — Modified NDWI

• But : mieux séparer l’eau des zones urbaines ou bâti.
• Formule :

.

.

• SWIR1 = bande courte longueur d’onde (B11)
• Valeurs : positives = eau ; améliore la détection par rapport au NDWI.
• Usage : suivi des plans d’eau en zone urbaine ou côtière.

BSI — Bare Soil Index

• But : détecter les sols nus et les zones dénudées.
• Formule :

.

.

• Blue = bande bleue (B02)
• Red = B04, NIR = B08, SWIR1 = B11
• Valeurs : plus BSI est élevé, plus la surface est dégagée.
• Usage : suivi des sols nus, zones urbaines, zones minérales.

NDBI — Normalized Difference Built-up Index

• But : détecter les zones bâties ou imperméabilisées.
• Formule :

.

.

• Valeurs : positives = bâti, négatives = végétation ou eau.
• Usage : cartographie urbaine, suivi de l’extension des villes.

EVI — Enhanced Vegetation Index

• But : améliorer la sensibilité de NDVI dans les zones très denses en végétation et corriger l’effet atmosphérique.
• Formule :

.

.

• Blue = B02, Red = B04, NIR = B08
• Valeurs : -1 à 1, similaire à NDVI mais plus robuste en zones denses.
• Usage : forêts, zones tropicales, suivi de la santé de la végétation.

SAVI — Soil Adjusted Vegetation Index

• But : mesurer la végétation en corrigeant l’effet du sol (utile dans les zones clairsemées).
• Formule :

.

.

où :

• NIR = bande proche infrarouge (Sentinel-2 B08)
• Red = bande rouge (B04)
• L = facteur de correction du sol, souvent 0.5
• Valeurs : -1 à 1, comme le NDVI, mais moins sensible au sol nu.
• Usage : végétation faible ou clairsemée (prairies, savanes, zones côtières avec faible couverture végétale).

UI — Urban Index (ou Urbanization Index)

• But : détecter et quantifier les zones urbaines.
• Formule simple :

.

.

• SWIR = bande courte longueur d’onde (Sentinel-2 B11 ou B12)
• NIR = proche infrarouge (B08)
• Valeurs : positives pour les zones imperméabilisées ou construites, faibles/négatives pour végétation ou eau.
• Usage : suivi de l’urbanisation, détection de bâtiments ou surfaces artificielles.

RDI — Redness Difference Index

• But : détecter les sols nus ou les zones riches en oxyde de fer (« rouge ») dans le paysage.
• Formule :

.

.

• Red = bande rouge (B04), Green = bande verte (B03)
• Valeurs : plus le sol est « rouge », plus le RDI est élevé.
• Usage : suivi de l’érosion, sols nus, zones minérales ou arides.

Tableau récapitulatif complet

Palette d’indices Sentinel

Chaque indice est présenté avec :

• son nom et acronyme,
• type de surface ciblé (végétation, eau, sol, urbain),
• valeurs typiques / échelle,
• couleur représentative pour visualisation sur carte (par exemple NDVI → vert, NDWI → bleu, BSI → marron clair, etc.).

Pourquoi utiliser un masque nuages ?

Les nuages et leurs ombres posent problème en imagerie satellite car :

1. Ils obscurcissent les surfaces : sol, végétation, eau, coraux.
2. Ils perturbent les indices (NDVI, NDWI…) car les valeurs radiométriques ne correspondent pas à la surface réelle.
3. Les zones couvertes de nuages peuvent générer des valeurs aberrantes : très hautes ou très basses, rendant les calculs de tendances faussés.

Donc, pour l’analyse d’indices et la détection de caractéristiques au sol ou dans l’eau, il est important de masquer ces pixels.

Masque nuages pour Sentinel-2 à 10 m

Sentinel-2 propose plusieurs options pour détecter les nuages à sa résolution native (10 m, 20 m ou 60 m selon la bande) :

a) Masque officiel fourni par ESA

• Fichier QA60 dans les produits Sentinel-2 L2A : masque binaire 10 m indiquant nuages et cirrus.
• Avantages : fiable et simple à appliquer.
• Inconvénients : masque grossier → ne différencie pas toujours bien nuages fins ou ombres.

b) Méthodes basées sur les indices spectrales

• Bande SWIR et NIR : les nuages sont très réfléchissants dans le visible mais moins dans le SWIR.
• Exemple d’index simple :

.

.

• Pixels clairs dans le visible et faibles dans SWIR → probablement des nuages.
• Permet de créer un masque nuages adapté au contexte, utile si on veut détecter nuages et cirrus locaux.

c) Algorithmes automatiques

• Fmask ou s2cloudless : détection sophistiquée utilisant machine learning et spectres multispectraux.
• Pratique pour grandes séries temporelles.

Comment on a procédé dans le script Colab

Dans lescript, on n’a pas utilisé le QA60 officiel, mais plutôt une détection par indices spectral et texture locale, adaptée aux zones côtières et à la super-résolution 1 m :

Étapes principales :

1. Normalisation des bandes :
   blue_n = blue / max
green_n = green / max
red_n = red / max
nir_n = nir / max
swir1_n = swir1 / max → toutes les valeurs entre 0 et 1 pour homogénéité.
2. Condition “nuage initial” basée sur 4 critères :
   • Très clair dans le visible (albedo_vis > 0.35)SWIR1 élevé (swir1_n > 0.15)Ratio bleu/rouge élevé (blue_n / red_n > 1.2)NIR élevé (nir_n > 0.25)

   → pixels respectant au moins 3/4 critères sont considérés comme nuages.

3. Masque “foam” pour les zones marines :

   • On détecte la mousse / déferlement des vagues via l’écart-type local sur la bande bleue et l’albédo visible.
   • On retire ces pixels du masque nuages :
     cloud_final = cloud_init & (~foam_mask)
   • Cela évite de masquer la zone corallienne ou la mer agitée.

4. Filtrage des petites composantes :

   • Supprime les nuages très petits qui pourraient être bruit → on garde uniquement les objets > 500 pixels.

5. Résultat final :

   • cloud_mask = masque nuages nettoyé, prêt à être utilisé pour exclure les pixels nuageux lors du calcul des indices.

Pourquoi cette méthode ?

• La scène est super-résolue (1 m) → le QA60 officiel n’existe pas à cette résolution.
• La zone est côtière / corallienne → mousse et déferlement doivent être distingués des nuages.
• L’approche spectral + texture locale permet :

  • détecter les nuages
  • ne pas masquer l’eau ou la mousse
  • travailler sur des images super-résolues générées par S2DR3

Mode opératoire

Ouvrez Google Colab: https://colab.research.google.com

Dans Collab, liez votre notebook à votre Google Drive

Maintenant, pour que le traitement ne prenne pas trop de temps, vous devez activer un GPU dans l’environnement d’exécution. Regardez dans le coin inférieur droit de votre notebook, vous verrez quelque chose comme « Python 3 » avec une icône de puce. Cliquez dessus et sélectionnez Changer le type d’environnement d’exécution. Une fenêtre s’ouvrira dans laquelle vous devrez rechercher Accélérateur matériel et sélectionner GPU T4. Enregistrez les modifications et le notebook redémarrera automatiquement avec la nouvelle configuration.

Copiez ce code dans une nouvelle cellule de votre bloc-notes et exécutez-le :

script Colab

# =========================================================
#   S2DR3 + Calcul indices Sentinel-2 + Masque nuages
#   Étude des coraux à Blue Bay / Île aux Aigrettes
# =========================================================

# --- 1. Monter Google Drive ---
from google.colab import drive
drive.mount('/content/drive')

# Crée un dossier pour les résultats
!mkdir -p /content/drive/MyDrive/Sentinel2_Coraux_S2DR3/output
!ln -s /content/drive/MyDrive/Sentinel2_Coraux_S2DR3/output /content/output

# --- 2. Installer le paquet S2DR3 ---
!pip -q install https://storage.googleapis.com/0x7ff601307fa5/s2dr3-20250905.1-cp312-cp312-linux_x86_64.whl

# --- 3. Importer le module principal ---
import s2dr3.inferutils
import os, glob, numpy as np, rasterio, scipy.ndimage as ndi
from skimage.transform import resize

# --- 4. Définir la zone d'intérêt et la date ---
lonlat = (57.73, -20.44)
date = '2025-10-09'

# --- 5. Lancer le traitement S2DR3 ---
s2dr3.inferutils.test(lonlat, date)

# =========================================================
#   6. Lecture du fichier _MS.tif
# =========================================================
root_dir = '/content/output'
files = sorted(glob.glob(os.path.join(root_dir, '**/*_MS.tif'), recursive=True),
               key=os.path.getmtime, reverse=True)
if not files:
    raise FileNotFoundError("❌ Aucun fichier *_MS.tif trouvé dans /content/output")
ms_path = files[0]
outdir = os.path.dirname(ms_path)
print(f" Fichier trouvé : {os.path.basename(ms_path)}")
print(f" Dossier : {outdir}")

# Lecture des bandes
with rasterio.open(ms_path) as src:
    profile = src.profile
    n_bands = src.count
    print(f"Nombre de bandes détectées : {n_bands}")
    bands = [src.read(i).astype('float32') for i in range(1, n_bands+1)]

def get_band(index, name):
    try:
        return bands[index-1]
    except IndexError:
        print(f" Bande {name} manquante, remplacée par zéro.")
        return np.zeros_like(bands[0])

# Attribution bandes principales Sentinel-2
blue  = get_band(1, "B02")
green = get_band(2, "B03")
red   = get_band(3, "B04")
nir   = get_band(4, "B08")
swir1 = get_band(9, "B11")
swir2 = get_band(10, "B12")

# =========================================================
#   7. Calcul des indices
# =========================================================
def safe_div(a, b, mask=None):
    res = np.where((b!=0) & (~np.isnan(a)) & (~np.isnan(b)), a/b, 0)
    if mask is not None:
        res[mask] = np.nan
    return res

L = 0.5  # facteur SAVI

indices = {
    'NDVI' : safe_div(nir - red, nir + red),
    'NDWI' : safe_div(green - nir, green + nir),
    'MNDWI': safe_div(green - swir1, green + swir1),
    'BSI'  : safe_div((swir1 + red) - (nir + blue), (swir1 + red) + (nir + blue)),
    'NDBI' : safe_div(swir1 - nir, swir1 + nir),
    'EVI'  : 2.5 * safe_div((nir - red), (nir + 6*red - 7.5*blue + 1)),
    'SAVI' : safe_div((nir - red) * (1 + L), (nir + red + L)),
    'UI'   : safe_div(swir1 - nir, swir1 + nir),
    'RDI'  : safe_div(red - green, red + green)
}

# =========================================================
#   8. Masque nuages léger adapté 1 m
# =========================================================
eps = 1e-8
blue_n, green_n, red_n, nir_n, swir1_n = [x/(np.nanmax(x)+eps) for x in [blue, green, red, nir, swir1]]
albedo_vis = (blue_n + green_n + red_n)/3

cond_bright = albedo_vis > 0.35
cond_swir = swir1_n > 0.15
cond_ratio = (blue_n/(red_n+1e-6)) > 1.2
cond_nir = nir_n > 0.25
cloud_init = (cond_bright.astype(int) + cond_swir.astype(int) + cond_ratio.astype(int) + cond_nir.astype(int)) >= 3

MNDWI_calc = np.where((green + swir1)!=0, (green - swir1)/(green + swir1), -1)
water_mask = MNDWI_calc > 0

size = 7
mean = ndi.uniform_filter(blue_n, size=size)
mean_sq = ndi.uniform_filter(blue_n**2, size=size)
local_std = np.sqrt(np.maximum(0, mean_sq - mean**2))
foam_mask = water_mask & (albedo_vis>0.25) & (local_std>0.03)

cloud_mask = cloud_init & (~foam_mask)
label_im, nb = ndi.label(cloud_mask)
sizes = ndi.sum(cloud_mask, label_im, range(1, nb+1))
mask_keep = np.zeros(nb+1, dtype=bool)
mask_keep[1:] = sizes >= 500
cloud_mask = mask_keep[label_im].astype(bool)

print(f" Nuages finaux : {cloud_mask.sum()} pixels")

# =========================================================
#   9. Sauvegarde indices individuels
# =========================================================
for name, data in indices.items():
    out_path = os.path.join(outdir, f"{name}.tif")
    profile_idx = profile.copy()
    profile_idx.update(dtype='float32', count=1, compress='lzw', nodata=None)
    with rasterio.open(out_path, 'w', **profile_idx) as dst:
        dst.write(data.astype('float32'), 1)
    print(f" {name} enregistré → {out_path}")

# =========================================================
#   10. Création du stack multibande robuste
# =========================================================
all_layers = list(indices.values()) + [cloud_mask.astype('float32')]
band_names = list(indices.keys()) + ['CLOUD_MASK']

ref_shape = indices['NDVI'].shape
for i, layer in enumerate(all_layers):
    if layer.shape != ref_shape:
        print(f" Redimensionnement couche {band_names[i]} de {layer.shape} → {ref_shape}")
        all_layers[i] = resize(layer, ref_shape, order=0, preserve_range=True, anti_aliasing=False).astype('float32')

profile_stack = profile.copy()
profile_stack.update(
    count=len(all_layers),
    dtype='float32',
    compress='lzw',
    nodata=None,
    height=ref_shape[0],
    width=ref_shape[1]
)

stack_path = os.path.join(outdir, 'indices_stack.tif')
with rasterio.open(stack_path, 'w', **profile_stack) as dst:
    for i, layer in enumerate(all_layers, start=1):
        dst.write(layer.astype('float32'), i)
        dst.set_band_description(i, band_names[i-1])

if os.path.exists(stack_path):
    print(f"\n Fichier multibande créé avec succès : {stack_path}")
    print(" Bandes : " + ', '.join(band_names))
else:
    print(f"\n Échec de création de {stack_path}")

print("\n Traitement complet terminé.")

Une fois exécuté vous aurez les fichiers .tif dans votre google drive

Conclusion : vers un traitement intégré des images Sentinel-2

L’intégration du calcul des indices spectraux directement dans le processus de super-résolution représente une évolution naturelle du traitement des données Sentinel-2.
En automatisant à la fois la génération de l’image à 1 m et la production des indices, on obtient un jeu de données homogène, complet et immédiatement exploitable pour l’analyse spatiale.

Appliquée aux milieux côtiers de Maurice et de Rodrigues, cette méthode ouvre de nouvelles perspectives :

• suivi de la végétation littorale et des mangroves,
• observation des platiers récifaux et des herbiers marins,
• détection des zones de turbidité ou d’envasement,
• cartographie des changements d’occupation du sol dans les zones côtières sensibles.

En réunissant toutes ces étapes dans un seul flux reproductible, ce script devient un outil de recherche et de suivi environnemental particulièrement adapté aux îles tropicales.
Il démontre aussi que les outils libres (Sentinel-2, Google Colab, QGIS) permettent aujourd’hui d’atteindre un niveau d’analyse spatiale autrefois réservé aux plateformes professionnelles.

Prochaine étape : intégrer ces indices dans une analyse chronologique multi-dates, afin de suivre l’évolution saisonnière ou post-cyclonique des habitats côtiers.

Si cet article vous a intéressé et que vous pensez qu'il pourrait bénéficier à d'autres personnes, n'hésitez pas à le partager sur vos réseaux sociaux en utilisant les boutons ci-dessous. Votre partage est apprécié !

Atilio Francois