Processing gdalogr:cliprasterbyextent en PyQGIS para cortar un ráster por extensión

En esta ocasión se va a probar en PyQGIS el método de processing gdalogr:cliprasterbyextent para cortar un ráster por extensión. Para el "corte" se va a emplear la extensión de un vectorial tipo polígono; tal como se presenta en la imagen siguiente:

Los puntos de la imagen anterior hacen referencia, respectivamente, a los puntos (xmin, ymin) y (xmax, ymax) correspondientes a la extensión del vectorial y cuya visualización ha sido posible mediante el plugin QuickWKT de QGIS.

Los parámetros que requiere la ejecución del método, 13 en total, se exponen en las dos imágenes siguientes.

Debido al gran número de parámetros a especificar, el código a continuación está debidamente documentado para facilitar su comprensión.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

import processing registry = QgsMapLayerRegistry.instance() vector_clip = registry.mapLayersByName('vector_clip') raster_output = registry.mapLayersByName('raster_output') extent = vector_clip[0].extent() xmin, ymin, xmax, ymax = extent.toRectF().getCoords() path = '/home/zeito/pyqgis_data/new_raster_clipped.tif' type = raster_output[0].type() provider = raster_output[0].dataProvider() dataType = provider.dataType(type) QGIS_TYPE = {QGis.Byte:'Byte', 2:'UInt16', 3:'Int16', 4:'UInt32', 5:'Int32', 6:'Float32', 7:'Float64'} new_type = QGIS_TYPE[dataType] RTYPE = { 'Byte': 0, 'Int16': 1, 'UInt16': 2, 'UInt32':3, 'Int32':4, 'Float32':5, 'Float64':6 } processing.runalg('gdalogr:cliprasterbyextent', raster_output[0], #INPUT <ParameterRaster> '', #NO_DATA <ParameterString> "%f,%f,%f,%f" % (xmin, xmax, ymin, ymax), #PROJWIN <ParameterExtent> RTYPE[new_type], #RTYPE <ParameterSelection> ---> 5 - Float32 4, #COMPRESS <ParameterSelection> ---> 4 - DEFLATE 75.0, #JPEGCOMPRESSION <ParameterNumber> 6.0, #ZLEVEL <ParameterNumber> 1.0, #PREDICTOR <ParameterNumber> False, #TILED <ParameterBoolean> 0, #BIGTIFF <ParameterSelection> ---> 0- False, #TFW <ParameterBoolean> None, #EXTRA <ParameterString> path) #OUTPUT <OutputRaster>

Ejecutado el código anterior en la Python Console de QGIS, al cargar la capa ráster "cortada" por la extensión del vectorial, se observa en la imagen siguiente que le faltan las dos últimas filas. No obstante, el ráster resultante está perfectamente alineado con el ráster base.

Este bug aparente puede ser resuelto calculando los índices de (fila, columna) del primer y último píxel de corte para establecer de manera más precisa los valores máximos y mínimos del ráster resultante.