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).
Best Answer
I managed to find a workaround for this issue by converting the data format to Int16 from Float32. The minimum value is then -32768 and can be processed as a nodata value. The following command did the trick:
There's probably a better solution, but this solves my immediate problem at least.