• Ir a navegación principal
  • Ir al contenido principal
  • Ir a la barra lateral primaria
  • Ir al pie de página
Logo web Geoinnova

Territorio Geoinnova - SIG y Medio Ambiente

El Blog de SIG, Territorio y Medio Ambiente

  • Asociación
    • Quienes somos
    • Geolibrería
    • Ofertas de empleo
    • Alianzas estratégicas
  • Consultoría
    • Consultoría y Desarrollo en Sistemas de Información Geográfica
    • Consultoría de Medio Ambiente
    • Turismo sostenible
    • Planificación urbana y ordenación del territorio
    • Planes Urbanos de Actuación Municipal
    • Geomarketing
    • QElectricGIS
  • Formación
    • Master GIS
    • Todos los cursos
  • Educación Ambiental
    • Organización de Campamentos
    • Organización de Itinerarios Ambientales
  • Coworking
  • Blog
    • SIG
    • Medio Ambiente
    • Teledetección
    • Programación y Desarrollo
    • Corporativo
  • Contacto
  • Inicio Blog
  • SIG
  • Programación y Desarrollo SIG
  • Teledetección
  • Medio Ambiente
  • Corporativo
Teledetección

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

Feb 24, 2021 · Por RemOT Technologies 3 comentarios

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, en este caso sobre el cálculo de tasseled cap con imágenes Sentinel-2.

La transformación Tasseled Cap es una combinación lineal de información espectral en un conjunto menor de bandas, normalmente tres ejes. A través de estas tres bandas se simplifican los componentes espectrales y permite tratarlas como variables físicas: brillo, verdor y humedad.

¿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:

BRI <- 0.3510*B02+ 0.3813*B03
          + 0.3437*B04+ 0.7196*B08
          + 0.2396*B11 + 0.1949*B12
  #Tasselled Cap - brightness

GRE <- (-0.3599* B02) + (-0.3533*B03)
          +(-0.4734*B04) + 0.6633*B08 
          +0.0087*B11 + (-0.2856*B12)
  #Tasselled Cap – greenness

WET <- (0.2578*B02 + 0.2305*B03 
          +0.0883*B04 + 0.1071*B08 
          +(-0.7611*B11) + (-0.5308*B12))
  #Tasselled Cap - wetness

Los coeficientes Tasselled Cap utilizados han sido obtenidos para datos de reflectancia Sentinel-2 y combinación de seis bandas (azul, verde, rojo, NIR, SWIR1 y SWIR2). Los coeficientes obtenidos han sido previamente comparados y validados con datos de Landsat-8 y MODIS (Shi y Xu, 2019).

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:

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

Bibliografía:

  • Index DataBase
  • 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
  • Shi, T., & Xu, H. (2019). Derivation of tasseled cap transformation coefficients for Sentinel-2 MSI at-sensor reflectance data. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 12(10), 4038-4048.

Este artículo ha sido actualizado a fecha 24 de Febrero de 2021

Artículos relacionados

Interfaz principal de CloudCompare. Nube de puntos LiDAR
10 motivos para aprender CloudCompare si trabajas con LiDAR
Mar 24, 2021
Teledetección
LiDAR
potree
Los 7 mejores visores LiDAR gratuitos
Feb 18, 2021
Teledetección
LiDAR
Imagen del SIOSE de Alta Resolución. Fuente: Geoportal SIOSE.
El nuevo SIOSE de Alta Resolución
Feb 16, 2021
Medio Ambiente, Teledetección
RemOT Technologies

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.

Interacciones con los lectores

Comentarios

  1. AvatarRicardo Alvarez P. dice

    Nov 19, 2020 en 14:10

    Gracias! Interesante útil.

    Responder
  2. AvatarSergio Morell dice

    Ene 20, 2020 en 17:10

    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.

    Responder
    • RemOT TechnologiesRemOT Technologies dice

      Ene 20, 2020 en 17:50

      Hola Sergio,
      Tienes razón, hay un error en la redacción del artículo y al copiar los valores se pasó eliminar esos dígitos. Procedemos a corregirlo.
      En cuanto a la segunda cuestión, decidimos seguir el planteamiento de https://www.indexdatabase.de ya que el algoritmo está calculado para Sentinel-2, y decidimos no alterar la fórmula dada. No obstante, es una buena cuestión y en este artículo utilizan la banda 11 (https://www.sciencedirect.com/science/article/pii/S0303243418309462).
      Muchas gracias por tu aporte.

      Responder

Deja una respuesta Cancelar la respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

Barra lateral primaria

Suscríbete al boletin

    Si te apetece puedes dejarnos tu nombre para que los correos tengan un trato personalizado

    Mediante el envío de mis datos personales confirmo que he leído y acepto la política de privacidad

    RESPONSABLE: Asociación Geoinnova. FINALIDAD: envío de publicaciones, promociones e información sobre cursos y eventos. LEGITIMACIÓN: tu legítimo consentimiento. DESTINATARIOS: Active Campaign con titular Active Campaign LLC, alojada en EEUU y suscrita al EU PrivacyShield. DERECHOS: acceso, rectificación, limitación, supresión de los datos (en [email protected])y a presentar reclamación ante una autoridad de control. INFORMACIÓN ADICIONAL: Política de privacidad.

    La último

    Python para ArcGIS

    Usando Python para ArcGIS: ArcPy

    Los Objetivos de Desarrollo Sostenible están orientados a conseguir un mundo sostenible en todos los aspectos.

    ¿Qué son los ODS u Objetivos de Desarrollo Sostenible?

    Ejemplo de clase de la API PyQGIS de QGIS

    PyQGIS para la programación en QGIS con Python

    huella de carbono

    Huella de carbono aplicada a planes y proyectos sometidos a evaluación ambiental

    Clonar repositorio Laravel desde Github

    Cómo clonar un proyecto de Laravel desde GitHub

    Cursos SIG y MA Asociación Geoinnova

    Agenda de cursos para el mes de Abril 2021

    Interfaz principal de CloudCompare. Nube de puntos LiDAR

    10 motivos para aprender CloudCompare si trabajas con LiDAR

    Footer

    Sobre Nosotros

    Territorio Geoinnova pretende ofrecer noticias y formación sobre el medio ambiente y las tecnologías de la información geográfica con interés a dichos sectores profesionales.

    Legal

    • Nota Legal

    El blog y todo su contenido se encuentra bajo una licencia de Creative Commons Reconocimiento-NoComercial-CompartirIgual 3.0 Unported.

    Todas las categorías

    • Actualidad Ambiental (86)
    • Alimentación Saludable (17)
    • Corporativo (124)
    • Formación (81)
    • Geolibrería (40)
    • Medio Ambiente (957)
    • Patrocinado (13)
    • Programación y Desarrollo SIG (22)
    • SIG (495)
    • Tecnología (69)
    • Teledetección (25)

    © 2021 · Desarrollada por Juan María Arenas - OikosMSP