Cálculo de Tasseled Cap para imágenes Sentinel-2 con R

Una semana más volvemos con un nuevo ejemplo sobre las múltiples posibilidades que nos ofrece R para la realización de análisis espacial automatizado.

¿Qué son los índices espectrales?

El uso de índices de teledetección está ampliamente extendido para el análisis de la superficie.

De un tiempo a esta parte, se ha hecho muy común, por ejemplo, el uso del Índice de vegetación de diferencia normalizada, más conocido como NDVI, aunque se viene utilizando desde el año 1973. Este índice es el más extendido y pone en relación los canales infrarrojo cercano (IRC) y rojo (R) para darnos información sobre el vigor de la vegetación, independientemente de la cantidad de cosecha o la calidad como en algunas ocasiones se vende.

Pero no sólo se utiliza el NDVI para el análisis de la vegetación o de la superficie.

El índice NDVI se basa en operaciones algebraicas realizadas sobre sólo dos bandas espectrales, la correspondiente a R e IRC. Sin embargo, Sentinel-2 es capaz de detectar información en 13 bandas, que también contienen información de interés en el estudio de las cubiertas vegetales y que permiten obtener más información.

Tasseled Cap

Kauth, R. J. and Thomas, G. S. (1976) pensaron en una transformación de la información contenida en las bandas que generase nuevos ejes que maximizasen la información contenida para utilizarla, inicialmente, en el campo de la agricultura. Se refirieron a dicha transformación como “Tasseled Cap” (TC).

Se trata de un método de compresión para reducir múltiples datos espectrales, concretamente de 6 bandas, en unas pocas bandas que permite comprender fenómenos importantes del desarrollo de cultivos en el espacio espectral, obteniendo los siguientes neocanales de información sobre:

• Brightness (brillo): está asociado a las variaciones de reflectancia del suelo
• Greeness (verdor): está correlacionado con el vigor de la vegetación
• Wetness (Humedad): está influído por las bandas en el IR medio y tiene que ver con la humedad vegetal y del suelo.

Limitaciones de Tasseled Cap

Una de las mayores limitaciones de la transformación Tasseled-cap es que solo se puede aplicar a sensores para los que ya se han desarrollado coeficientes de transformación. Actualmente, esto incluye Landsat, MODIS, ASTER, Ikonos, SPOT y, además, Sentinel 2, tal y como veremos a continuación.

SCRIPT PARA CALCULAR TASSELED CAP EN IMÁGENES SENTINEL-2

Las imágenes Sentinel-2A y Sentinel-2B se pueden descargar en formato JPEG2000 con extensión .jp2. Se trata de una compresión utilizada para imágenes grandes y de alta calidad.

Vamos a crear una lista de todas las rutas que apuntan a las bandas disponibles de nuestra imagen a través de la función list.files(). Indicaremos nuestra ruta a la carpeta descargada con los datos Sentinel como path.

path <- "D://MSIL1C_20191202T105421_N0208_R051_T30TXM_20191202T112708"

Vamos a definir una búsqueda a través del patrón de caracteres “jp2”, a través del parámetro pattern para que solamente se seleccionen los archivos con esa extensión. Vamos a preservar la ruta completa al archivo a través del parámetro full.names = TRUE y vamos a indicar que realice una búsqueda recursiva en las subcarpetas de nuestra ruta inicial con recursive = TRUE.

jp2_files <- list.files(path = path, pattern = ".jp2", full.names = TRUE, recursive = TRUE)

El objeto jp2_files va a contener las bandas de nuestra imagen. Sobre este objeto, que solamente contiene la ruta a cada una de las bandas, vamos a utilizar la función raster() para leer la imagen en formato raster. Para ello, debemos cargar previamente la librería raster:

library(raster)

Vamos a ayudarnos de la función subset() y grepl() para seleccionar la ruta a la banda específica que en este caso contiene la cadena de texto “B02” y que corresponde con la banda del Azul. Y sobre esta selección aplicaremos la función raster().

B02 <- raster(subset(jp2_files, grepl("B02",jp2_files)))

Podemos previsualizar nuestra banda B02:

plot(B02)

Vamos a repetir este procedimiento con cada una de las bandas que vayamos a utilizar de ahora en adelante.

Para calcular nuestros Tasselled Cap, vamos a utilizar distintas bandas, en algunos casos, con distinta resolución espacial. En ese caso debemos utilizar la función aggregate() o disaggregate() para convertir nuestras bandas a una resolución mayor o menor, pero siempre común para poder operar nuestra matriz de datos celda por celda. Una vez obtenida la resolución común para todas las capas vamos a calcular nuestro Tasselled Cap con las siguientes fórmulas:

SBI <- 0.3037*B02+0.2793*B03
  +0.4743*B04+0.5585*B08+
  0.5082*B10+0.1863*B12
  #Tasselled Cap - brightness
  
GVI <- (-0.2848*B02)+(-0.2435*B03)+(-0.5436*B04)
  +0.7243*B08+0.0840*B11+(-0.1800*B12)
  #Tasselled Cap - vegetation
  
WET <- 0.1509*B02+0.1973*B03+0.3279*B04
  +0.3406*B08+(-0.7112*B11)+(-0.4572*B12)
  #Tasselled Cap - wetness

Vamos a recordar de nuevo que estas ecuaciones son específicas para Sentinel-2.

Finalmente, vamos a crear un objeto ráster multibanda con la función stack() en el que vamos a incluir las tres bandas de nuestro Tasselled Cap. Aplicaremos la función plot() para devolver el resultado de nuestra composición multibanda Tasselled Cap (brightness – greeness – wetness).

Esta es sólo una pequeña muestra de la multitud de aplicaciones que tiene R para el trabajo con imágenes satélite y Sistemas de Información Geográfica. No olvides echar un vistazo a los cursos que impartimos en colaboración con nuestros amigos de RemOT Technologies que profundizan en el uso de R:

• Curso de Análisis espacial con R
• Curso de Análisis de regresión y modelado espacial en R

Bibliografía:

• Kauth, R. J. and Thomas, G. S. (1976) – The tasselled cap – a graphic description of the spectraltemporal development of agricultural crops as seen by Landsat
• Index DataBase