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

1 estrella2 estrellas3 estrellas4 estrellas5 estrellas (8 votos, media: 5,00 sobre 5)
Cargando…

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)
Azul_con_reflectancia_escalados
Banda B02 – Azul con valores de reflectancia escalados

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

Composicion_RGB_Tasselled_Cap
Composición RGB del 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:

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
Artículo anteriorEl PDTI de Cullera elaborado por Geoinnova logra el máximo nivel de calificación
Artículo siguienteLa revolución de los geodatos
RemOT Technologies
RemOT Technologies es una empresa colaboradora de Geoinnova, especialista en desarrollo de aplicaciones GIS para la web, cartografía y análisis espaciales automatizados. Desarrolla visores web cartográficos, aplicaciones basadas en geolocalización y plugins y personalizaciones de QGIS para resolver problemas de distintos ámbitos, como por ejemplo la gestión de parcelas o la gestión de redes. RemOT ya ha trabajado en varios proyectos de desarrollo de la mano de Geoinnova y, además, también ha participado, entre otros proyectos como en el desarrollo del Atlas Nacional de España. Ha sido reconocida como una de las 100 mejores empresas geoespaciales del mundo en el año 2019 por Geoawesomeness, una de las 50 Startup españolas de futuro por la revista Emprendedores y posee el sello Pyme Innovadora del Ministerio de Ciencia Innovación y Universidades. En la parte dedicada a formación, Lucía Martínez, Marcos Gimeno y Miquel Febrer, imparten diversos cursos en la plataforma de Geoinnova.

2 Comentarios

  1. Hola,

    Creo que hay un error en el código.

    Al calcular las componentes SBI GVI y WET se utilizan coeficientes con 5 decimales. Parece que el último decimal sobra. Posiblemente este error se deba a la fuente de donde se extrajo la información: https://www.indexdatabase.de/db/i-single.php?id=91 El último número corresponde a la banda de Sentinel-2 que se debe utilizar, no es un dígito propio del coeficiente.

    Por otro lado, planteo otra cuestión.
    La banda 10 de Sentinel-2 (cirrus) aporta información atmosférica, no de la superficie terrestre. En esta parte del espectro (1375 nm) la atmósfera es opaca. Es por ello que al realizar la corrección atmosférica utilizando Sen2Cor para obtener la reflectancia de la superficie terrestre, el algorimo elimina la B10. ¿Sigue siendo útil entonces para calcular la componente SBI?

    Saludos,

    S.M.

Dejar respuesta

Please enter your comment!
Please enter your name here