Blog SIG & Territorios

En dos artículos anteriores, exploramos cómo pasar de imágenes Sentinel-2 estándar (resolución de 10 m) a una versión de alta resolución de 1 m utilizando el módulo S2DR3, desarrollado para la investigación sobre los corales de Mauricio. Este paso permite obtener un nivel de detalle sin precedentes a partir de datos abiertos, revelando estructuras finas que a menudo son invisibles a 10 m: praderas marinas, plataformas, zonas arenosas o incluso transiciones sutiles entre la vegetación y el arrecife.

Tradicionalmente, el cálculo de los índices espectrales (NDVI, NDWI, MNDWI, BSI, NDBI, etc.) se realiza en QGIS, una vez descargadas y procesadas las imágenes. Sin embargo, a menudo es más eficaz generar estos índices directamente en la fuente, es decir, desde la producción de las imágenes Sentinel-2 de alta resolución.

Este enfoque presenta varias ventajas:

• Automatización completa: los índices se calculan desde la creación del conjunto de datos, sin necesidad de procesamiento manual adicional en QGIS.
• Coherencia espacial: los índices utilizan exactamente las mismas correcciones radiométricas y la misma superresolución que las bandas originales.
• Ahorro de tiempo: sin exportaciones intermedias ni configuración manual de algoritmos ráster.
• Preparación para el análisis: los resultados se pueden utilizar inmediatamente en QGIS o en cualquier otra herramienta SIG.

El script que se presenta a continuación automatiza todo el flujo:

1. descarga y superresolución de la escena Sentinel-2,
2. generación de una máscara de nubes inteligente (con exclusión de las zonas de oleaje),
3. cálculo directo de los principales índices espectrales,
4. creación de un archivo multibanda listo para el análisis.

Esta integración simplifica considerablemente el paso de la imagen bruta a los mapas de interpretación, al tiempo que garantiza la reproducibilidad del tratamiento científico.

Los índices calculados

NDVI — Índice de vegetación de diferencia normalizada.

• Objetivo: medir la densidad y la salud de la vegetación.
• Fórmula:

.

.

• NIR = banda del infrarrojo cercano (B08)
• Red = banda roja (B04)
• Valores: -1 a 1; >0,2 vegetación, <0,1 suelo desnudo o agua.
• Uso: seguimiento de la vegetación, estimación de la biomasa, agricultura.

NDWI — Índice de agua normalizado por diferencia

• Objetivo: detectar la presencia de agua en el paisaje.
• Fórmula:

.

.

• Verde = banda verde (B03)
• NIR = banda infrarroja cercana (B08)
• Valores: positivos = agua, negativos = suelo/vegetación.
• Uso: cartografía de lagos, ríos, humedales.

MNDWI — NDWI modificado

• Objetivo: separar mejor el agua de las zonas urbanas o construidas.
• Fórmula:

.

.

• SWIR1 = banda de longitud de onda corta (B11)
• Valores: positivos = agua; mejora la detección en comparación con el NDWI.
• Uso: seguimiento de masas de agua en zonas urbanas o costeras.

BSI — Índice de suelo desnudo

• Objetivo: detectar suelos desnudos y zonas desnudas.
• Fórmula:

.

.

• Azul = banda azul (B02)
• Red = B04, NIR = B08, SWIR1 = B11
• Valores: cuanto más alto es el BSI, más despejada está la superficie.
• Uso: seguimiento de suelos desnudos, zonas urbanas, zonas minerales.

NDBI — Índice de diferencia normalizada de zonas construidas

Objetivo: detectar zonas construidas o impermeabilizadas.

Fórmula:

.

.

• Valores: positivos = construido, negativos = vegetación o agua.
• Uso: cartografía urbana, seguimiento de la expansión de las ciudades.

EVI — Índice de vegetación mejorado

• Objetivo: mejorar la sensibilidad del NDVI en zonas con vegetación muy densa y corregir el efecto atmosférico.
• Fórmula:

.

.

• Azul = B02, Rojo = B04, NIR = B08
• Valores: -1 a 1, similar al NDVI pero más robusto en zonas densas.
• Uso: bosques, zonas tropicales, seguimiento de la salud de la vegetación.

SAVI — Índice de vegetación ajustado al suelo

Objetivo: medir la vegetación corrigiendo el efecto del suelo (útil en zonas poco densas).

Fórmula:

.

.

donde:

• NIR = banda del infrarrojo cercano (Sentinel-2 B08)
• Red = banda roja (B04)
• L = factor de corrección del suelo, a menudo 0,5
• Valores: -1 a 1, como el NDVI, pero menos sensible al suelo desnudo.
• Uso: vegetación escasa o dispersa (praderas, sabanas, zonas costeras con poca cobertura vegetal).

UI — Índice urbano (o índice de urbanización)

Objetivo: detectar y cuantificar las zonas urbanas.

Fórmula simple:

.

.

• SWIR = banda de longitud de onda corta (Sentinel-2 B11 o B12)
• NIR = infrarrojo cercano (B08)
• Valores: positivos para zonas impermeables o construidas, bajos/negativos para vegetación o agua.
• Uso: seguimiento de la urbanización, detección de edificios o superficies artificiales.

RDI — Índice de diferencia de rojez

• Objetivo: detectar suelos desnudos o zonas ricas en óxido de hierro («rojo») en el paisaje.
• Fórmula:

.

.

• Rojo = banda roja (B04), Verde = banda verde (B03)
• Valores: cuanto más «rojo» es el suelo, mayor es el RDI.
• Uso: seguimiento de la erosión, suelos desnudos, zonas minerales o áridas.

Tabla resumen completa

Gama de índices Sentinel

Cada índice se presenta con:

• su nombre y acrónimo,
• tipo de superficie objetivo (vegetación, agua, suelo, urbano),
• valores típicos/escala,
• color representativo para su visualización en el mapa (por ejemplo, NDVI → verde, NDWI → azul, BSI → marrón claro, etc.).

¿Por qué utilizar una máscara de nubes?

Las nubes y sus sombras plantean un problema en las imágenes satelitales porque:

1. Oscurecen las superficies: suelo, vegetación, agua, corales.
2. Alteran los índices (NDVI, NDWI, etc.) porque los valores radiométricos no se corresponden con la superficie real.
3. Las zonas cubiertas por nubes pueden generar valores atípicos: muy altos o muy bajos, lo que distorsiona los cálculos de tendencias.

Por lo tanto, para el análisis de índices y la detección de características en el suelo o en el agua, es importante ocultar estos píxeles.

Máscara de nubes para Sentinel-2 a 10 m.

Sentinel-2 ofrece varias opciones para detectar nubes en su resolución nativa (10 m, 20 m o 60 m según la banda):

a) Máscara oficial proporcionada por la ESA.

• Archivo QA60 en los productos Sentinel-2 L2A: máscara binaria de 10 m que indica nubes y cirros.
• Ventajas: fiable y fácil de aplicar.
• Inconvenientes: máscara tosca → no siempre distingue bien entre nubes finas y sombras.

b) Métodos basados en índices espectrales

• Banda SWIR y NIR: las nubes son muy reflectantes en el visible, pero menos en el SWIR.
• Ejemplo de índice simple:

.

.

• Píxeles claros en el visible y débiles en SWIR → probablemente nubes.
• Permite crear una máscara de nubes adaptada al contexto, útil si se desea detectar nubes y cirros locales.

c) Algoritmos automáticos

• Fmask o s2cloudless: detección sofisticada que utiliza aprendizaje automático y espectros multiespectrales.
• Práctico para grandes series temporales.

Cómo se procedió en el script Colab

En el script no se utilizó el QA60 oficial, sino una detección por índices espectrales y textura local, adaptada a las zonas costeras y a la superresolución de 1 m:

Etapas principales:

1. Normalización de las bandas :
   blue_n = blue / max
green_n = green / max
red_n = red / max
nir_n = nir / max
swir1_n = swir1 / max → todos los valores entre 0 y 1 para homogeneidad.
2. Condición «nube inicial» basada en 4 criterios :
   • Muy claro en el visible (albedo_vis > 0,35) SWIR1 alto (swir1_n > 0,15) Relación azul/rojo elevada (blue_n / red_n > 1,2) NIR elevado (nir_n > 0,25)

   → los píxeles que cumplen al menos 3/4 criterios se consideran nubes.

3. Máscara «espuma» para las zonas marinas :

   • Se detecta la espuma/el rompimiento de las olas mediante la desviación estándar local en la banda azul y el albedo visible.
   • Se eliminan estos píxeles de la máscara de nubes :
     cloud_final = cloud_init & (~foam_mask)
   • Esto evita enmascarar la zona coralina o el mar agitado.

4. Filtrado de pequeños componentes :

   • Elimina las nubes muy pequeñas que podrían ser ruido → solo se conservan los objetos > 500 píxeles.

5. Resultado final :

   • cloud_mask = máscara de nubes limpia, lista para ser utilizada para excluir los píxeles nubosos al calcular los índices.

¿Por qué este método?

• La escena tiene una resolución muy alta (1 m) → el QA60 oficial no existe con esta resolución.
• La zona es costera/coralina → hay que distinguir el musgo y las olas de las nubes.
• El enfoque espectral + textura local permite:

  • detectar las nubes
  • no ocultar el agua ni el espuma
  • trabajar con imágenes de alta resolución generadas por S2DR3

Modo de funcionamiento

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

En Collab, vincule su cuaderno a su Google Drive.

Ahora, para que el procesamiento no lleve demasiado tiempo, debe activar una GPU en el entorno de ejecución. Mire en la esquina inferior derecha de su cuaderno, verá algo como «Python 3» con un icono de punto. Haga clic en él y seleccione Cambiar el tipo de entorno de ejecución. Se abrirá una ventana en la que deberá buscar Acelerador de hardware y seleccionar GPU T4. Guarde los cambios y el cuaderno se reiniciará automáticamente con la nueva configuración.

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é.")

Una vez ejecutado, tendrá los archivos .tif en su Google Drive.

Conclusión: hacia un procesamiento integrado de las imágenes Sentinel-2

La integración del cálculo de los índices espectrales directamente en el proceso de superresolución representa una evolución natural del procesamiento de los datos Sentinel-2.
Al automatizar tanto la generación de la imagen a 1 m como la producción de los índices, se obtiene un conjunto de datos homogéneo, completo y directamente utilizable para el análisis espacial.

Aplicado a los entornos costeros de Mauricio y Rodrigues, este método abre nuevas perspectivas:

• seguimiento de la vegetación litoral y los manglares,
• observación de los arrecifes y las praderas marinas,
• detección de zonas de turbidez o sedimentación,
• cartografía de los cambios en la ocupación del suelo en zonas costeras sensibles.

Al reunir todas estas etapas en un único flujo reproducible, este script se convierte en una herramienta de investigación y seguimiento medioambiental especialmente adecuada para las islas tropicales.
También demuestra que las herramientas libres (Sentinel-2, Google Colab, QGIS) permiten hoy en día alcanzar un nivel de análisis espacial que antes estaba reservado a las plataformas profesionales.

Próxima etapa: integrar estos índices en un análisis cronológico multidate, con el fin de seguir la evolución estacional o posciclónica de los hábitats costeros.

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