La temperatura de la superficie terrestre se puede estimar o calcular usando las bandas térmicas de Landsat 8. Simplemente se requiere aplicar un conjunto de ecuaciones a través de una calculadora de imágenes ráster (ArcMap, ArcGIS Pro, QGIS).
El primer paso consiste en descargar una imagen Landsat 8 de algún lugar en particular, descomprimir, y chequear cierta información necesaria (dentro de los metadatos) para ejecutar este procedimiento.
Este tutorial muestra cómo calcular la temperatura de superficie terrestre (LST, por sus siglas en inglés Land Surface Temperature) usando las bandas térmicas de Landsat 8. Particularmente, la banda 10 como banda térmica, y las bandas 4 y 5 para calcular el índice de vegetación de diferencia normalizada (NVDI, por sus siglas en inglés Normal Difference Vegetation Index).
Para calcular el LST, se usa las fórmulas del USGS (mayor información en el artículo Algorithm for Automated Mapping of Land Surface Temperature Using LANDSAT 8 Satellite Data), este ejemplo realiza simplemente los cálculos sin entrar en detalle sobre los fundamentos. A continuación se sintetiza el proceso en seis pasos:
1.- Cálculo de TOA (Top of Atmospheric) radiancia espectral.
TOA (L?) = ML * Qcal + AL
Dónde:
ML = factor de reescalamiento multiplicativo específico de banda (valor disponible en el archivo de metadatos MTL, en la línea RADIANCE_MULT_BAND_x, donde x es el número de banda).
Qcal = corresponde a la banda 10.
AL = factor de reescalamiento aditivo específico de la banda (valor disponible en el archivo de metadatos MTL, en la línea RADIANCE_ADD_BAND_x, donde x es el número de banda)
TOA = 0.0003342 * “Band 10” + 0.1
Por consiguiente se debe resolver la ecuación usando la herramienta Raster Calculator en ArcMap.
2.- Conversión de TOA a Brightness Temperature (Temperatura de brillo)
BT = (K2 / (ln (K1 / L?) + 1)) − 273.15
dónde:
K2 y K1 = constantes de conversión térmica específicas de la banda a partir de los metadatos.
L? = TOA
Por lo tanto para obtener los resultados en Celsius, la temperatura radiante se ajusta sumando el cero absoluto (aprox. -273,15°C).
BT = (1321.0789 / Ln ((774.8853 / “%TOA%”) + 1)) – 273.15
3.- Calcular el NDVI
NDVI = (Banda 5 – Banda 4) / (Banda 5 + Banda 4)
Ten en cuenta que el cálculo del NDVI es importante porque, posteriormente, se debe calcular la proporción de vegetación (Pv), que está altamente relacionada con el NDVI, y la emisividad (ε), que está relacionada con el Pv.
NDVI = Float(Banda 5 – Banda 4) / Float(Banda 5 + Banda 4)
4.- Calcular la proporción de vegetación Pv
Pv = Square ((NDVI – NDVImín) / (NDVImáx – NDVImín))
Por lo general los valores mínimos y máximos de la imagen NDVI se pueden visualizar directamente en la imagen (tanto en ArcGIS, QGIS, ENVI, Erdas Imagine), caso contrario se debe abrir las propiedades del ráster para obtener aquellos valores.
Pv = Square((«NDVI» – 0.216901) / (0.632267 – 0.216901))
5.- Calcular la Emisividad ε
ε = 0.004 * Pv + 0.986
Simplemente aplicar la fórmula en la calculadora ráster, el valor de 0.986 corresponde a un valor de corrección de la ecuación.
6.- Calcular la temperatura de superficie de la tierra
LST = (BT / (1 + (0.00115 * BT / 1.4388) * Ln(ε)))
Finalmente aplicar la ecuación de LST para obtener el mapa de temperatura de superficie.
Como resultado del proceso desarrollado se cuenta con un mapa de temperatura de superficie de la tierra, cabe señalar que no es igual a la temperatura del ambiente.
Muy buen contenido como siempre! Muchas gracias!
Una pregunta: para este caso es recomendable trabajar con las imágenes LANDSAT 8 ya con corrección atmosférica y radiométrica? o es mejor trabajarlas sin esta?
Estimado Andrés, claro lo puedes hacer tranquilamente.
Saludos,
Franz
Gracias por el aporte.
Una consulta, has tenido oportunidad de evaluar el LST plugin en QGIS?
El mismo funciona para L5, L7 y L8, opera de la siguiente manera:
Convierte DN a radiancia.
Luego radiancia a temperatura de brillo.
Calcula NDVI y luego estima la emisividad en superficie.
La estimacion de LST se realiza mediante aplicacion de ec de Planck.
Gracias
No lo he probado, gracias por compartir, ya lo chequeo.
franzpc, thanks! And thanks for sharing your great posts every week!
Buenas.
Me gustaría saber si una vez introducida la variable de Emisividades en la fórmula. ¿ estaríamos hablando de una temperatura REAL de superficie terrestre? o estamos hablando de temperatura APARENTE de superficie?
Muchas gracias
Son estimaciones, lo que obtendríamos es la temperatura aparente de la superficie terrestre.
¿Por qué solo utilizas la proporción de vegetación?
yo incluyo la proporción de suelo, de forma que la fórmula de la emisividad, teniendo en cuenta las reflexiones que pueda haber, queda así:
ε= εvegetación*Pvegetacion + εsuelo*(1-Pvegetacion)(1-1.74*Pvegetacion)+1.7372*Pvegetacion (1-Pvegetacion)
Gracias por tu aporte, lo voy a replicar, puesto que esta metodología es general, y no lo había visto de aquella forma.
Como calculo la Proporción del suelo?
Hola. Buen día Franzpc,
me gustaría saber si es posible utilizar la información generada como temperatura de un cultivo?
Personalmente no considero conveniente, puesto que estima la temperatura de la superficie terrestre.
El valor de 0.004 en la emisividad a que corresponde?
Buenos días tengo una duda, de donde se obtiene el valor 0.0003342 del TOA.
Muchas gracias.
Es un valor de radiancia (RADIANCE_MULT_BAND_10 RADIANCE_MULT_BAND_11) de las bandas 10 y 11 que se encuentra en los metadatos (ahí lo vas a encontrar como 3.3420E-04).
Para mayor información puedes chequear:
https://www.esri.com/arcgis-blog/products/product/analytics/deriving-temperature-from-landsat-8-thermal-bands-tirs/
http://www.un-spider.org/sites/default/files/LDCM-L8.R1.pdf
Buen día, para este procedimiento es necesario realizar algún tipo de corrección a la imágen?
… Otra consulta, tienes el procedimiento para imágenes) Landsat 7
Saludos
Buen dia Franz, tengo una duda, cuando intente sacar Pv funciona bien si la imagen esta recortada, pero si saco el NDVI a la imagen completa me genera el rango de NDVI de -1 a 1, cuando aplico la ecuacion para Pv en este caso, realiza la operacion pero no se ve ninguna imagen, esto se puede solucionar?
Si usas toda la imagen, probablemente tendrías que fijar como NoData los valores que limitan la imagen (fondo negro del recuadro).