Script de python para determinar sucesivamente NDVI, SAVI, LAI y εNB lanzado en cónsola de bash (Linux Debian)

Hace unos meses, establecí un procedimiento para determinar sucesivamente NDVI, SAVI, LAI y εNB mediante un script de bash empleando las bandas 3 y 4 de Landsat. Los valores de εNB (emisividad superficial) son usados en el cálculo de la temperatura superficial empleando la banda térmica (b6) de Landsat. En este caso, el script de bash se traduce y codifica al lenguaje de script de python, empleando las librerías de grass.script, el cual tiene la ventaja de que produce un código más legible y de ejecución sustancialmente más rápida que el de bash. No obstante, debido a las modificaciones del .bashrc y .grassrc6, es posible ejecutar el script de python desde cónsola de bash sin necesidad de tener abiertos GRASS o QGIS (cónsola de GRASS).

Para ello, se propuso el siguiente código que es válido para cualquier dataset que incluya las bandas 3 y 4 de Landsat renombradas como b3_corte y b4_corte (para diferenciarlas de las originales de mayor extensión). Es interesante notar que el uso de grass.mapcalc emplea el modificador $ antes del nombre del ráster y su referencia dentro de la instrucción si tiene que buscarlo dentro del dataset. Sin embargo, tal referencia es innecesaria si el ráster a ser usado ha sido previamente generado en tiempo de ejecución.

#!/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")

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

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

print "Hecho!"

print "Calculando eNBf ..."

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

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

print "Hecho!"

La ejecución del código produjo la siguiente salida en cónsola de bash:

cuyo tiempo es significativamente menor que el del script de bash. La imagen del ráster final de εNB fue la siguiente:

About these ads