Tratamiento de nubes en imágenes Sentinel con Python: estrategias efectivas

Los satélites de la serie Sentinel desarrollados por la Agencia Europea del Espacio (ESA) son un recurso muy valioso para el análisis geográfico de diversos aspectos relacionados con el medio ambiente, agricultura, desastres naturales… Sin embargo, la presencia de nubes en las imágenes plantean problemas. El más evidente, pero no único, es que bloquean la vista del suelo. En este artículo abordaremos varias estrategias para el tratamiento de nubes en imágenes Sentinel con Python.

El problema de las nubes en las imágenes de satélite

Las nubes bloquean la vista de la superficie terrestre, por lo que dificultan la observación de los elementos que deseamos monitorizar, como las masas de agua, la vegetación, o el uso del suelo. Pero además las nubes pueden afectar a la reflectancia medida por los satélites lo que influye en la calidad de algunas bandas. Esto tiene impacto si de lo que se trata es de obtener índices como el NDVI, porque puede distorsionar los resultados.

La medida más obvia es buscar imágenes de la zona de estudio que no tengan nubes. A este respecto es importante destacar que mediante el paquete PySTAC Client podemos acceder a catálogos STAC para obtener imágenes Sentinel con un umbral de cobertura de nubes definido.

Por ejemplo podemos buscar imágenes con cobertura menor del 10% en el catálogo.

La opción de trabajar con imágenes que no contengan nubes no siempre es posible, pues puede que necesitemos utilizar una imagen en una fecha o intervalo de tiempo concreto. En ese caso la limitación temporal puede implicar que no estén disponibles imágenes sin nubes.

1. Crear un máscara para eliminar las zonas cubiertas de nubes

Este es un método habitual que consiste en eliminar de la imagen aquellas zonas que contiene nubes. Podemos obtener las zonas cubiertas de nubes de varias formas.

1.1. Uso de la capa de clasificación de escena (SCL)

Empleo de SCL (Scene classification). La clasificación de escenas se desarrolló para distinguir entre píxeles nublados, píxeles claros y píxeles de agua de los datos de Sentinel-2 y es el resultado del algoritmo de clasificación de escenas de la ESA para productos de Level-2A. El mapa de clasificación se produce para cada producto Sentinel-2 Nivel 2A a una resolución de 20 y 60 m y los valores de bytes del mapa de clasificación se organizan como se muestra en la siguiente tabla:

Según esto las clases que utilizaremos para crear la máscara son la 9 y la 10, que como puedes observar son las que hacen referencia a las nubes. El código que vamos a utilizar utiliza NumPy y rasterio y es el siguiente:

import numpy as np
from rasterio.warp import reproject, Resampling

# Abrir la imagen RED
with rasterio.open("B04.tif") as src:
    red_band = src.read(1)
    profile = src.profile  # Guardar metadatos
    red_shape = red_band.shape  # (altura, anchura)

# Abrir la banda SCL (Scene Classification Layer)
with rasterio.open("SCL.tif") as scl_src:
    scl = scl_src.read(1)
    scl_profile = scl_src.profile  # Guardar metadatos de SCL

# Crear una máscara de nubes (clases 9 y 10 en SCL)
cloud_mask = (scl == 9) | (scl == 10)

# Reescalar la máscara SCL a la resolución de la banda RED
scl_resampled = np.zeros(red_shape, dtype=scl.dtype)  # Crear array vacío
reproject(
    source=scl,
    destination=scl_resampled,
    src_transform=scl_src.transform,
    dst_transform=src.transform,
    src_crs=scl_src.crs,
    dst_crs=src.crs,
    resampling=Resampling.nearest  # Reescalado con el vecino más cercano
)

# Aplicar la máscara a la banda RED

red_band = red_band.astype(np.float32)  # Convertir a float32
cloud_mask_resized = (scl_resampled == 9) | (scl_resampled == 10)  # Recalcular máscara
red_band[cloud_mask_resized] = np.nan  # Sustituir píxeles nublados por NaN

Lo primero es importar las librerías:

• numpy: Se usa para manejar matrices de datos (imágenes en formato raster).
• rasterio.warp: Permite reproyectar y reescalar imágenes con diferentes resoluciones.

Lo que vamos a hacer es eliminar las nubes de una banda, en este caso de la banda roja (B04) de Sentinel. Por lo tanto abrimos la imagen y obtenemos sus metadatos (profile) y las dimensiones de la banda roja, porque las vamos a necesitar cuando hagamos la re-proyección.

Luego abrimos la banda SCL (Scene Classification Layer).

Mediante; cloud_mask = (scl == 9) | (scl == 10) creamos una máscara booleana en la que los píxeles con valores 9 o 10 (nubes) se marcan como True, y el resto como False.

Reescalar las imágenes y aplicar la máscara

Es necesario reescalar porque las imágenes B04.tif (banda roja) y SCL.tif (clasificación de escenas) tienen diferentes resoluciones espaciales. La banda roja (B04) tiene una resolución de 10 metros por píxel, y la banda de clasificación SCL: 20 metros por píxel. Es decir, cada píxel en la SCL representa un área mayor que los píxeles de la banda roja. Si no se ajustara la resolución de la SCL a la de la banda roja, los tamaños de las matrices serían diferentes y no podríamos aplicar la máscara correctamente.

Una vez re-escalado podemos crear la máscara. Primero se convierte la banda roja a float32 para poder manejar valores NaN.

A continuación se emplea: cloud_mask_resized = (scl_resampled == 9) | (scl_resampled == 10)

que tiene como objetivo identificar los píxeles con nubes en la imagen re-escalada (scl_resampled), de manera que se puedan enmascarar en la banda roja (B04). Se utilizan los valores 9 y 10 que como indicamos entes se corresponden con las nubes. Por último se convierten las zonas con nubes a NaN.

Guardar la imagen corregida

Para finalizar guardamos la imagen que obtenido sin las nubes:

profile.update(dtype=rasterio.float32)
with rasterio.open("red_sin_nubes.tif", "w", **profile) as dst:
    dst.write(red_band.astype(np.float32), 1)

Visualizar el resultado

Para comparar la nueva imagen con la imagen de la banda roja inicial podemos emplear Matplotlib:

Para obtener esta representación hemos utilizado este código:

import numpy as np
import rasterio
import matplotlib.pyplot as plt

# Cargar la imagen original (Banda Roja)
with rasterio.open("B04.tif") as src:
    red_band = src.read(1)

# Cargar la imagen sin nubes
with rasterio.open("red_sin_nubes.tif") as src:
    red_sin_nubes = src.read(1)

# Crear una figura con 2 gráficos (1 fila, 2 columnas)
fig, axes = plt.subplots(1, 2, figsize=(12, 6))

# Mostrar la imagen original
ax1 = axes[0]
im1 = ax1.imshow(red_band, cmap="Reds")  # Colormap de tonos rojos
ax1.set_title("Imagen Original (Banda Roja)")
ax1.axis("off")  # Ocultar ejes
fig.colorbar(im1, ax=ax1, fraction=0.046, pad=0.04)  # Agregar barra de colores

# Mostrar la imagen sin nubes
ax2 = axes[1]
im2 = ax2.imshow(red_sin_nubes, cmap="Reds")
ax2.set_title("Imagen Sin Nubes (Banda Roja)")
ax2.axis("off")
fig.colorbar(im2, ax=ax2, fraction=0.046, pad=0.04)

# Mostrar la comparación
plt.tight_layout()
plt.show()

1.2. Detección de nubes con índices espectrales

Una alternativa al método anterior es el empleo de algoritmos para calcular índices espectrales que nos sirvan para identificar las nubes en las imágenes.

El índice de nieve de diferencia normalizada (NDSI) en Sentinel-2 se obtiene como la relación de dos bandas: una en el VIR (Banda 3) y otra en el SWIR (Banda 11). Aunque es un índice para obtener la nieve también puede ser de aplicación en la detección de nubes. Las nubes tienen alta reflectancia en las bandas azul y SWIR, pero baja en el rojo e infrarrojo cercano (NIR).

Los valores de NDSI obtenidos con la ecuación anterior cuando son superiores a 0,4 pueden indicar nubes.

Este método que puede ser efectivo para detectar nubes gruesas puede confundir otras nubes con nieve y agua. En la vista anterior vemos el resultado de aplicar el índice NDSI sobre la banda roja (B04) siendo el resultado poco satisfactorio, lo que implicaría hacer ajustes o modificar el umbral del valor de 0,4 que se a utilizado en este caso.

2. Tratamiento de las nubes mediante Machine Learning

Desde el punto de vista de Python que es el que estamos siguiendo aquí, hay varias herramientas que nos pueden ayudar a tratar el problema de las  nubes. Los siguientes métodos de Machine Learning son de aplicación:

• Random Forest o SVM entrenados con píxeles etiquetados como nubes/no-nubes.
• Redes neuronales (CNNs) para segmentación automática.

En cuanto a las librerías de Machine Learning que podemos utilizar para hacer estas operaciones tenemos:

• Scikit-learn es una biblioteca de aprendizaje automático (Machine Learning) para Python. Se utiliza para entrenar modelos de clasificación, regresión, clustering y reducción de dimensión. Esta librería nos proporciona los métodos Random forest o SVM. Mediante estas técnicas podemos entrenar un modelo que sea capaz de identificar las nubes.
• TensorFlow y Keras son bibliotecas de Deep Learning que permiten entrenar redes neuronales para tareas como clasificación, segmentación y detección de objetos.

Después de identificar las nubes, se puede interpolar los valores vecinos para rellenar los huecos. Sin embargo esta interpolación presenta riesgos. Al interpolar los valores estamos suponiendo que los «huecos a rellenar» tienen las mismas características que las áreas colindantes. Esta suposición no tiene por qué cumplirse siempre y es especialmente problemáticas en zonas de alta variabilidad donde podemos encontrarnos con un mosaico de superficies (zonas con agua, áreas urbanas, tierras de cultivo…)

3. Emplear imágenes de diferentes fechas

Otra estrategia que podemos emplear, es «rellenar» las zonas cubiertas con nubes con datos de otras imágenes de satélite de la misma zona pero que han sido obtenidas en una fecha libre de nubes. Esta estrategia tiene la ventaja de que se obtiene una imagen sin nubes y sin pérdida de información. Pero tiene el inconveniente de que se necesitan otras imágenes sin nubes lo que no siempre es posible.