lunes, 30 de octubre de 2017

Usando pygeoprocessing.geoprocessing.raster_calculator

El módulo pygeoprocessing de python tiene la opción de códificar una herramienta de cálculo ráster (raster_calculator) de una manera mucho más amigable que la de PyQGIS y equiparable a las de Processing. Como es "third part" puede instalarse con pip o easy_install.

Si importan el módulo en la Python Console (ver código siguiente), un help(geop.raster_calculator) les dará información detallada acerca de los parámetros a usar.

El código a continuación permite determinar el NDVI a partir de las bandas NIR y RED correspondientes. No hace falta tener cargadas en la Map View de QGIS las bandas respectivas.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
import pygeoprocessing.geoprocessing as geop
from osgeo import gdal

path = '/home/zeito/pyqgis_data/'

listraster_uri2 = [(path + 'b3_ref.tif',1), (path + 'b4_ref.tif',1)]

rasterout_uri2 = path + 'NDVI_geo.tif'

def ndvi(r1,r2):
   result_ndvi = (r2-r1)/(r1 + r2)
   return result_ndvi

geop.raster_calculator(base_raster_path_band_list = listraster_uri2,
                       local_op = ndvi, 
                       target_raster_path = rasterout_uri2,
                       datatype_target = gdal.GDT_Float32,
                       nodata_target = -999,
                       calc_raster_stats = False)

Después de la ejecución del código anterior en la Python Console de QGIS, el ráster resultante se cargó en la Map View y se colorizó con un archivo de estilo *qml que segrega las clases correspondientes a cuerpos de agua, suelo y vegetación. El resultado fue el siguiente:


No hay comentarios: