Suma de rásteres con pygeoprocessing evitando la generación de nodata values

Generalmente, cuando se utiliza cualquier versión de Raster Calculator (processing, pygeoprocessing, etc) proveniente del software libre uno de los requerimientos es que los rásteres tengan la misma resolución y estén perfectamente alineados.

Cuando se piensa en comenzar una rutina de trabajo con la rasterización de un vectorial como máscara, la limitación que se tiene es que, a pesar de tener control sobre la resolución, los rásteres resultantes no estarán alineados entre si para diferentes features y se producen nodata values que dificultan el álgebra de mapas en aquellas regiones donde aparecen.

Para evitar este inconveniente de la rasterización se puede usar en su lugar la opción de 'clip by mask' de Processing, para el ráster base, la cual garantiza alineación e idéntica resolución para todos los rásteres producidos.

El script personalizado (Processing) siguiente permite el corte por máscara produciendo para cada feature un ráster binario donde los valores son uno para la máscara y cero más allá de esta. Por tanto, en el área donde se solapan todos los rásteres la suma es máxima.

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

##Raster=group ##InputRaster=raster ##mask=vector ##Clipped=output raster mask = processing.getObject(mask) xmin = mask.extent().xMinimum() xmax = mask.extent().xMaximum() ymin = mask.extent().yMinimum() ymax = mask.extent().yMaximum() mask_extent = str(xmin)+ ',' + str(xmax)+ ',' +str(ymin)+ ',' +str(ymax) path = processing.runalg('gdalogr:cliprasterbymasklayer', InputRaster, #INPUT <ParameterRaster> mask, #MASK <ParameterVector> "", #NO_DATA <ParameterString> False, #ALPHA_BAND <ParameterBoolean> False, #CROP_TO_CUTLINE <ParameterBoolean> False, #KEEP_RESOLUTION <ParameterBoolean> 5, #RTYPE <ParameterSelection> 4, #COMPRESS <ParameterSelection> 75, #JPEGCOMPRESSION <ParameterNumber> 6, #ZLEVEL <ParameterNumber> 1, #PREDICTOR <ParameterNumber> False, #TILED <ParameterBoolean> 0, #BIGTIFF <ParameterSelection> False, #TFW <ParameterBoolean> None, #EXTRA <ParameterString> Clipped) #OUTPUT <OutputRaster>

En la imagen siguiente se tienen las tres máscaras vectoriales sobre el ráster base. Los solapamientos se observan porque se ha dado 50 % de transparencia a las dos primeras capas.

El resultado de utilizar el script anterior se encuentra en la imagen siguiente. También se observan los solapamientos porque, al igual que para los vectoriales, se ha asignado el mismo nivel de transparencia (50 %) a los dos primeros rásteres.

Finalmente, se usa el script siguiente que emplea el pygeoprocessing raster calculator para obtener la suma de rásteres.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

import pygeoprocessing.geoprocessing as geop from osgeo import gdal path = '/home/zeito/pyqgis_data/' listraster_uri = [(path + 'clipped_raster1.tif',1), (path + 'clipped_raster2.tif',1), (path + 'clipped_raster3.tif',1)] rasterout_uri= path + 'sum_clipped_rasters.tif' def sum(r1,r2,r3): result_sum=(r1+r2+r3) return result_sum geop.raster_calculator(base_raster_path_band_list=listraster_uri, local_op = sum, target_raster_path=rasterout_uri, datatype_target = gdal.GDT_Byte, nodata_target = None, calc_raster_stats = False)

En la imagen a continuación puede observarse el ráster resultante. Con la ayuda del plugin Value Tool de QGIS se ha corroborado que la suma es 3, tal como se esperaba, donde todos ellos se solapan (area blanca).