[GIS] QGIS Raster calculation with nodata values

Question

I am using QGIS 2.18.13.

I have a few shapefiles. From inside QGIS, using python's processing.runalg function and supplying "grass7:v.to.rast.attribute", I converted those to raster files which have value of 1 everywhere shapes were present and "no data" (255) beyond them.

Those rasters share the same spatial location and I would like to "add them" so if in the same location there are 3 values 1 in final raster there will be value 3. So, I will be happy to have nodata converted to 0s for that purpose.

I tried following things where I didn't get results, I found some limitations or I got no results at all.

1. Raster -> raster calculator – I don't have any control over
   resolution. It was produced uncompressed tiff (400MB! from my
   10000×10000 raster), 32 bit precision and nodata values changed
   to -3.40282e+38. Where it was nodata before I get nodata so,
   resulting raster is some kind of intersection rather than sum.

2. Toolbox -> GRASS GIS 7 -> Raster -> r.null. I get this error: "The following layers were not correctly generated. NullRaster"

3. toolbox -> GRASS GIS 7 -> Raster -> r.mapcalc. I don't know how to write an expression. I always get either parse error or:

   r.mapcalc expression="abc = 1" --overwrite 
   r was unexpected at this time.
   

   there are usually 4 lines in the log like this:

   ERROR: Sorry is not a valid option 
   

4. toolbox -> SAGA -> Raster calculus -> Raster calculator – that one seems not to throw errors, however, I am not sure how to handle NaN's?

I would like to add this to python script so I can automate the process if my base datasets change.

Answer 1

You are going to have issues if you try to follow your workaround. They pop up because whatever version of raster calculator that you can use, usually, it needs same resolution and aligning for each rasterized shapefile for an optimal performance. In this case, it's preferable to employ 'clip by mask' by using a base raster with all values equal 1.

For trying out my approach, I used three shapefiles and base raster of following image. You can observe each overlapping because it was assigned 50 % of transparency in two first layers.

For clipping by mask each shapefile I used following personalized processing script.

##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>

Result of running above script for each shapefile it can be visualized at following image; where it was also assigned 50 % of transparency in two first layers to observe each raster overlapping.

Finally, I used following python script, that uses pygeoprocessing raster calculator to get sum clipped rasters.

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)

At following image it can be observed resulting raster. With Value Tool QGIS plugin it was corroborated that sum was 3, as expected, where three rasters are overlapped (white area).