1 LAS MÁSCARAS

Las nubes y las condiciones atmosféricas constituyen un obstáculo cuando se trabaja con imágenes de satélite. Tanto nubes como sus sombras pueden dificultar el análisis de muchas escenas debido a que sus valores de reflectancia son muy altos (las nubes son muy brillantes, ya que dispersan toda la radiación que vuelve al sensor) o demasiado oscuros (las sombras representan radiación que ha sido bloqueada o absorbida). Aunque las nubes pueden enmascararse, interpolarse (sustituirse por datos de otra fecha) o convertirse en un tipo de cobertura separada del resto en una clasificación, como regla general es conveniente evitar el uso de escenas con mucha nubosidad.

En el caso de que las nubes ocupen una superficie reducida en la escena, una alternativa es la utilización de una máscara (capa) en la que se identifican los píxeles que probablemente correspondan a nubes o sombras. Una vez hecho esto, se pueden cambiar los valores de esos píxeles a NA en las bandas que se utilizarán para posteriores análisis.

Metadatos de la máscara de nubes raster

Existen varios procedimientos para generar esta máscara de nubes. En todos los casos deben seguirse una recomendaciones generales comunes:

Este capítulo está basado fundamentalmente en la siguiente página web

1.1 Generación de la máscara de nubes con RStoolbox

library(raster)
## Loading required package: sp
library(RStoolbox)
library(ggplot2)

Importamos el ejemplo y vemos que tiene nubes en la zona oriental

data(lsat) 
ggRGB(lsat, 
      stretch = "lin")

La obtención de la máscara de nubes requiere varios pasos.

  • PASO 1. Cálculo de un índice de nubosidad. Este índice se calcula a partir de las bandas azul y TIR. Una propiedad que distingue a las nubes es una radiación térmica baja (las nubes son frías por que están emplazadas a gran altitud) y su alta reflectancia en la región visible del espectro. Es esta información la que puede utilizarse para identificar y enmascarar automáticamente las nubes, si bien hay que tener cuidado porque la nieve puede tener propiedades radiativas similares.
cldmsk  <- cloudMask(lsat, 
                     blue = 1, 
                     tir = 6)
ggR(cldmsk, 
    2, 
    geom_raster = TRUE)

En segundo lugar, definimos un umbral, reutilizando el índice calculado previamente. . Además, se calcula alrededor de los píxeles con nubes una función de crecimiento.

cldmsk_final <- cloudMask(cldmsk, 
                          threshold = 0.1,       # Todos los valores superiores a ese umbral serán enmascarados
                          buffer = 5)            # Número de píxeles usados como buffer. 

ggRGB(lsat, stretch = "lin") +
   ggR(cldmsk_final[[1]], 
       ggLayer = TRUE, forceCat = TRUE, geom_raster = TRUE) +
   scale_fill_manual(values = c("red"), na.value = NA)
## Warning: Removed 88592 rows containing missing values (geom_raster).

Estimación del desplazamiento de la sombra de las nubes

Se puede realizar interactivamente o no interacticamente. Interactivamente: click en el pixel de la nube y en el/los píxeles correspondientes a las sombras. ## Not run: shadow <- cloudShadowMask(lsat, cldmsk_final, nc = 2)

No interactivamente: debe añadirse una estimación del movimiento de la sombra (shiftEstimate) en unidades del mapa

shadow <- cloudShadowMask(lsat,                                 # Imagen original
                          cldmsk_final,                         # Máscara de nubes creada previamente
                          shiftEstimate = c(-16,-6))            # Estimación del movimiento de la sombra 

csmask <- raster::merge(cldmsk_final[[1]], 
                        shadow)

ggRGB(lsat, stretch = "lin") +
        ggR(csmask, ggLayer = TRUE, forceCat = TRUE, geom_raster = TRUE) +
        scale_fill_manual(values = c("blue", "yellow"), 
        labels = c("shadow", "cloud"), na.value = NA)
## Warning: Removed 88214 rows containing missing values (geom_raster).

1.2 Generación de la máscara de nubes con ficheros de Landsat.

Los archivos comprimidos con imágenes de Landsat incorporan a menudo ficheros que brindan información adicional sobre el error y la incertidumbre en los datos. Entre ellos aparecen ficheros que continen una máscara de nubes, en la que se identifican qué píxeles corresponderán a nubes o sombras generadas por estas nubes. En el caso de Landsat 8, aparecen 2 ficheros con la etiqueta cfmask.

  • cfmask

  • cfmask_conf

El archivo cfmask.tif contiene varios números. Para conocer su significado es necesario abrir el fichero de metadatos (data/week-08/landsat/LC80340322016189-SC20170128091153/LC80340322016189LGN00.xml) desde el explorador de windows.

Metadatos de la máscara de nubes raster

De acuerdo con este fichero, el valor 2 corresponde a las sombras de las nubes, mientras que el 4 corresponde a las nubes.

Para cargar las imágenes necesarias para esta actividad, además de establecer el directorio de trabajo con la función setwd, la función list.files lista todos los archivos de la carpeta deseada que contengan la combinación de letras band y que acaben con la extensión .tif(el signo del dólar al final).

setwd("D:/G174_OCW/Lab_3_preprocesamiento_imagenes_satelite/")
ejemplo_nubes <- list.files("./datos/mascara_nubes/ficheros_nubes/",
                            pattern = glob2rx("*band*.tif$"), 
                            full.names = TRUE)

Con los ficheros de la lista se creará primero una RasterStack y, a partir de él, un RasterBrick.

ejemplo_nubes_st <- stack(ejemplo_nubes)
ejemplo_nubes_br <- brick(ejemplo_nubes_st)

A continuación representamos gráficamente esa imagen

par(col.axis = "white", col.lab = "white", tck = 0)               # Color del eje blanco, sin marcas
plotRGB(ejemplo_nubes_br,
        r = 4, g = 3, b = 2,
        stretch = "hist",
        main = "Imagen color natural con nubes",
        axes = TRUE)
box(col = "white")                                                # Cambia el cuadro a blanco para que no haya borde en su trama

A continuación, también se puede representar la capa de máscara que corresponde a las nubes.

mascara_nubes_conf <- raster("./datos/mascara_nubes/ficheros_nubes/LC80340322016189LGN00_cfmask_conf_crop.tif")

plot(mascara_nubes_conf, 
     main = "Máscara de nubes")

A continuación, también se puede representar la capa de máscara que corresponde a las sombras. El valor 2 corresponde a las sombras, mientras que el 4 corresponde a las nubes. El 1 corresponde a superficies acuáticas.

mascara_nubes <- raster("./datos/mascara_nubes/ficheros_nubes/LC80340322016189LGN00_cfmask_crop.tif")

plot(mascara_nubes,
  main = "Máscara de nubes y sombras")

1.2.1 Crear capa de máscara

En este caso, desea establecer todos los valores mayores que 0 en la máscara de ráster en NA.

par(xpd = FALSE, 
    mar = c(0, 0, 1, 5))

Creación de una máscara que combina nubes y sombras. Todos los valores mayores que 2 (coresponden a nubes y sombras) son transformados en NA

mascara_nubes[mascara_nubes > 1] <- NA

Representar gráficamente la máscara

plot(mascara_nubes,
     main = "Máscara de nubes y sombras",
     col = c("green"),
     legend = FALSE,
     axes = FALSE,
     box = FALSE)
par(xpd = TRUE)                                     # Se obliga a la leyeda a aparecer fuera de la imagen
legend(x = mascara_nubes@extent@xmax, 
       y = mascara_nubes@extent@ymax,
       c("Not masked", "Masked"),
       fill = c("green", "white"),
       bty = "n")

Todos los píxeles de color VERDE corresponden a píxeles no enmascarados, mientras que los píxeles en blanco están enmascarados por que corresponden a nubes y sombras.

1.2.2 Aplicar la máscara

Una vez obtenida la máscara de nubes, se puede superponer al RasterBrick para enmascarar todos los píxeles que corresponden a nubres.

ejemplo_nubes_mascara <- mask(ejemplo_nubes_br,            # RasterBrick que queremos "limpiar" de nubes
                               mask = mascara_nubes)               # Fichero de máscara creado anteriormente

Comprobación de si el fichero con la máscara está en memoria del ordenador

inMemory(ejemplo_nubes_mascara)
## [1] TRUE
class(ejemplo_nubes_mascara)
## [1] "RasterBrick"
## attr(,"package")
## [1] "raster"

Representar gráficamente la imagen enmascarada en color natural

# Primero se cambian todos los ejes al color blanco y se desactivan las etiquetas
par(col.axis = "white", col.lab = "white", tck = 0)

# Dibujo de los datos
plotRGB(ejemplo_nubes_mascara,
        r = 4, g = 3, b = 2,
        main = "Imagen en color natural sin estiramiento",
        axes = TRUE)
box(col = "white")

Observe anteriormente que no tuve que usar la función de estiramiento para forzar que los datos se trazaran en R. Esto se debe a que los píxeles extremadamente brillantes que representaban nubes ahora se eliminan de sus datos.

# Primero cambie todos los ejes al color blanco y desactive las marcas
par(col.axis = "white", col.lab = "white", tck = 0)

# Dibujar los datos
plotRGB(ejemplo_nubes_mascara,
        r = 4, g = 3, b = 2,
        stretch = "lin",
        main = "Imagen en color natural con estiramiento",
        axes = TRUE)
box(col = "white")

Eliminamos todos los objetos emplazados en el Global Environment

rm(list = ls())

1.3 Reemplazar los valores de los píxeles de la máscara de nubes por otros datos

Una posibilidad para el análisis de imágenes de satélite parcialmente cubierts por nubes consiste en reemplazar los píxeles cubiertos por nubes (en blanco) por los valores de otro ráster no afectado por este problema. Puede ser una escena tomada días antes o días después. El procedimiento también requiere varios pasos.

  • PASO 1: la escena de Landsat “limpia” debe convertirse en un RasterBrick.
no_nubes <- list.files("./datos/mascara_nubes/ficheros_limpios/",
           pattern = glob2rx("*band*.tif$"),                        
           full.names = TRUE)   
  
no_nubes_st <- stack(no_nubes)
no_nubes_br <- brick(no_nubes_st)

Representación gráfica

plotRGB(no_nubes_br,
        r = 4, g = 3, b = 2,
        stretch = 'lin')

Ahora importaremos las imágenes con nubes. La siguiente sintaxis realiza una lista con todos los archivos que contengan la palabra banda y su extensión sea .tif.

nubes <- list.files("./datos/mascara_nubes/ficheros_nubes",
                    pattern = glob2rx("*band*.tif$"),
                    full.names = TRUE)

Creación de un raster raster stack y un raster brick con los ficheros de la lista

nubes_st <- stack(nubes)
nubes_br <- brick(nubes_st)

Representación gráfica del raster brick

plotRGB(nubes_br,
        r = 4, g = 3, b = 2,
        stretch = "lin")

Aplicación de la máscara de nubes a los datos con nubes.

cloud_mask <- raster("./datos/mascara_nubes/ficheros_nubes/LC80340322016189LGN00_cfmask_crop.tif")
cloud_mask[cloud_mask > 1] <- NA

nubes_br_mask <- mask(nubes_br,
                      mask = cloud_mask)

plotRGB(nubes_br_mask,
        r = 4, g = 3, b = 2,
        stretch = "lin")

Verifique si los rásteres están en la misma extensión y CRS.

compareRaster(no_nubes_br, 
              nubes_br_mask)

Las extensiones son diferentes, recortemos una a la otra.

no_nubes_br <- crop(no_nubes_br, extent(nubes_br))

compareRaster(no_nubes_br, 
              nubes_br_mask)
## [1] TRUE

La función raster::cover() reemplaza cada píxel con valor NA en el fichero con nubes por el valor de reflectancia del mismo píxel y banda procedentes del raster sin nubres.

raster_limpio <- cover(nubes_br_mask,
                       no_nubes_br)

plotRGB(raster_limpio, 
        r = 4, g = 3, b = 2, 
        stretch = 'lin')

1.4 Exportar el fichero raster

Como es frecuente cuando se realizan diferentes tareas de preprocesamiento, cada una de las bandas puede adquirir nombres complicados.

raster_limpio
## class      : RasterBrick 
## dimensions : 177, 246, 43542, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 455655, 463035, 4423155, 4428465  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : LC8034032//band1_crop, LC8034032//band2_crop, LC8034032//band3_crop, LC8034032//band4_crop, LC8034032//band5_crop, LC8034032//band6_crop, LC8034032//band7_crop 
## min values :                     0,                     0,                     0,                     0,                     0,                     0,                     0 
## max values :                  2735,                  3001,                  3923,                  4219,                  4985,                  4524,                  4491

Por eso, es conveniente que las bandas adopten un nombre más sencillo.

names(raster_limpio) <- paste("B",
                    1:7,
                    sep="")
raster_limpio
## class      : RasterBrick 
## dimensions : 177, 246, 43542, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 455655, 463035, 4423155, 4428465  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      :   B1,   B2,   B3,   B4,   B5,   B6,   B7 
## min values :    0,    0,    0,    0,    0,    0,    0 
## max values : 2735, 3001, 3923, 4219, 4985, 4524, 4491

Finalmente, el objeto ráster es exportado con formato Geotiff.

writeRaster(raster_limpio, 
            filename = "./datos/mascara_nubes/LC8034032_sin_nubes.tif",
            bylayer=TRUE,                                             # Crea un fichero separado para cada capa
            format="GTiff", 
            overwrite=TRUE)