Script de Python (consola de GRASS-QGIS) para determinar sucesivamente NDVI, SAVI, LAI, εNB y temperatura superficial (Ts) en Ubuntu

Hace unos meses establecí un procedimiento para determinar sucesivamente NDVI, SAVI, LAI y εNB mediante un script de python empleando las bandas 3 y 4 de Landsat 5. Los valores de εNB (emisividad superficial) son usados en el cálculo de la temperatura superficial (Ts) empleando la banda térmica (b6) de Landsat. Según el manual del SEBAL se determinan mediante la siguiente expresión (ecuación de Plank modificada):

temp

donde Rc es la radiancia térmica corregida obtenida a partir de la radiancia térmica de la superficie y K1 (1260.56 ºK) y K2 (607.76 Wm-2sr-1μm-1; las mismas unidades de la radianza) son constantes. Como el cálculo de Rc requiere de información de difícil obtención se va a prescindir de ésta en la determinación de Ts usando en su lugar los valores de radianza (L). El manual del Sebal advierte que el uso de radianzas no corregidas podría producir temperaturas subestimadas hasta en 5 grados para las zonas más calientes de la imagen.

El script propuesto es el siguiente:

#!/usr/bin/env python
# -*- coding: utf-8
 
import grass.script as grass
from os import system
 
system("clear")
 
print "Calculando NDVI ..."
 
grass.mapcalc("NDVI=($b4_corte-$b3_corte)/($b4_corte+$b3_corte+0.)",
                b3_corte="b3_corte",
                b4_corte="b4_corte",
                overwrite="true",
				quiet="true")
 
print "Hecho!"
 
print "Calculando SAVI con L=0.1 ..."
 
grass.mapcalc("SAVI_0.1= (($b4_corte-$b3_corte)/($b4_corte+$b3_corte+0.1))*(1+0.1)",
                b3_corte="b3_corte",
                b4_corte="b4_corte",
                overwrite="true",
				quiet="true")

print "Hecho!"
 
print "Calculando LAI..."
 
grass.mapcalc("LAI=if(SAVI_0.1 < 0.1, 0.00001,(if(0.1 < SAVI_0.1 && SAVI_0.1 < 0.687,-log((0.69-SAVI_0.1)/0.59)/0.91,if(SAVI_0.1 > 0.687,6,0))))",
                overwrite="true",
				quiet="true")

print "Hecho!"
 
print "Calculando eNB ..."
 
grass.mapcalc("eNB=if(LAI<3 && NDVI>0.,0.97+0.0033*LAI,if(LAI>=3 && NDVI>0.,0.98))",
                overwrite="true",
				quiet="true")

print "Hecho!"

print "Calculando eNBf ..."
 
grass.mapcalc("eNBf=if(eNB == 0, 0.99, eNB)",
                overwrite="true",
				quiet="true")

print "Hecho!"

print "Calculando Ts ..."

grass.mapcalc("Ts = 1260.56/log(eNBf*607.76/b6.rad2+1)",
				overwrite="true",
				quiet="true")

print "Hecho!"

Su ejecución produjo la siguiente batería de imágenes:

ndvi_savi

lai_enb

thermal_Ts

Esta entrada fue publicada en Código Python, GRASS, QGIS, SIG, Software Libre. Guarda el enlace permanente.

9 respuestas a Script de Python (consola de GRASS-QGIS) para determinar sucesivamente NDVI, SAVI, LAI, εNB y temperatura superficial (Ts) en Ubuntu

  1. Adhemar Conde dijo:

    Hola Jose,

    Estoy intentando utilizar el scrip de Python que has realizado para el calculo de la Temperatura superficial, pero tengo un problema al ejecutar el script, la consola me arroja un error:

    import grass.script as grass
    ImportError: No module named grass.script

    Hasta donde he podido averiguar (espero no estar equivocado) “grass.script” es un conjunto de librerías de Grass para Python, pero no sé como conseguir esta librería. No sé si esta librería está incluida en el paquete de instalación de Grass Gis o en el paquete de instalación de Qgis.

    En fin, si tengo instalado en mi ordenador Grass Gis o Qgis, el script que tu has elaborado debería correr sin ningún problema? O es que debo instalar esta librería como se instala la librería de GDAL?

    Tal vez este tema es un tema sin mayores complicaciones, y la pregunta que realizo es muy básica, pero debo comentarte que sé muy poco de este tema, soy nuevo utilizando estas herramientas.

    Muchas gracias por una futura respuesta.

    Saludos cordiales,
    Adhemar.

    • Supongo que te estás refiriendo a Ubuntu. Si ese error pasa verifica (con synaptic) si tienes instalado qgis-plugin-grass. Ese script lo implemente en Ubuntu 12.04 Wubi emulado en un Windows 7.

  2. Adhemar Conde dijo:

    Hola Jose,

    Gracias por tu respuesta. He seguido tu consejo y si, si tengo instalado qgis-plugin-grass.

    Quisiera pedirte un favor, me enseñarías a utilizar la consola de Grass-Qgis? Estoy trabajando sobre Ubuntu 14.04 y tengo Qgis 2.2.

    1. En Qgis habilito el complemento GRASS
    2. Presiono el botón “Abrir directorio de mapas”, aquí ya tengo creado la “Base de datos (Gisdbase)”, la “Localizacion” y el “directorio de mapas”. En el “directorio de mapas” tengo dos imágenes que corresponden a la banda 3 (Red) y la banda 4 (NIR).
    3. Después de haber cargado estos directorios se me habilitan los demás iconos, presiono en el icono “abrir herramientas de GRASS”.
    4. Después utilizo la opción “shell-Consola de GRASS”
    5. El script que has elaborado lo he copiado en un archivo .py y lo guardo dentro de una carpeta (con el nombre GDAL) dentro de Documentos.
    6. En la consola utilizo el comando: cd Documentos/GDAL/
    7. Una vez dentro de la carpeta GDAL (en la consola Grass-Qgis) ejecuto el comando:
    python Ts.py (Ts.py es el script que tu has elaborado)

    Después de todo esto, y con los nervios que todo salga bien, me arroja un error:

    Calculando NDVI …
    ERROR: G_getenv(): La variable LOCATION_NAME no se ha establecido.
    ERROR: G_getenv(): La variable LOCATION_NAME no se ha establecido.

    No sé como establecer LOCATION_NAME. No sé como ejecutar el script, no sé donde deberían estas las imágenes.

    Siento haberme extendido. Y nuevamente muchas gracias por la ayuda prestada.

    Saludos,
    Adhemar.

  3. Adhemar dijo:

    Hola José,

    Un par de preguntas por favor.

    1. Al calcular la Ts utilizas la banda térmica del landsat (b6). Esta banda esta corregida, es decir, has determinado las temperaturas (°K) en el tope de la atmósfera. Es esta banda 6 corregida la que utilizas para el cálculo de la Ts, no es cierto?

    2. Cual será el valor mínimo y máximo de la Ts que has calculado con la imagen satelital?. Me imagino que sus valores deben ser parecidos a los valores de la banda 6 corregida.

    Muchas gracias.
    Saludos,
    Adhemar.

    • 1. No. Es con los valores de radianza corregida. Date cuenta que es la misma fórmula de Plack. Si usas las radianzas no corregidas entonces los valores son algo menores (hasta 5 grados en las zonas más calientes).

      2. Los valores, para esta región de Perú, oscilaron entre 8 y 30 grados centígrados (suma 273,15 para grados kelvin). Los valores de Ts serían algo mayores que la llamada temperatura de brillo. Como yo no corregí las radianzas entonces esta Ts que determiné es la misma temperatura de brillo.

  4. José, porqué es difícil la obtención de información para Rc.

    Gracias..

  5. Le agradecería pueda ayudarme con las siguientes inquietudes.

    Estoy elaborando una metodología para poder determinar la Temperatura Superficial con imágenes Landsat 8.
    Lo que no e podido entender es para Landsat 8, que tiene 2 bandas térmicas (10 y 11). como saber con cual de las bandas trabajar.

    Lo mismo sucede con el sensor ETM+ donde también viene 2 bandas, una de alta y baja ganancia, como explicar el uso de una de estas bandas y el porque.

    Le estaré muy agradecido desde ya.

    Saludos.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s