2

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 10000x10000 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.

  • No data is not the same that NaN (Not a Number). Normally, no data is a number used for this purpose (e.g., -999). NaN is obtained as result of a not allowed operation [a/0, log(0), sqrt(-a), a + NaN, etc] or considering values "out of extent", e.g., for first pixel in a raster five pixels are considered NaN for its eight possible neighbors ("border effect"). – xunilk Oct 28 '17 at 3:37
  • I appreciate the difference. Therefore it is even more so annoying that there is lack of consistency. In general tab under layer properties, it is stated that there is No-Data Value of 255. In style when applying gradient it tries to apply from 'nan' to 1, even though I know layer contains only 1s and 255s. when I load all values it correctly sets minimum and maximum at 1. So if there is a value (of 255) why can't it be altered? When I run raster calculator (let's say raster - 55) In return I get -54s, and No-Data (of -3.40282e+38). So adding two rasters will work on intersection only... – Bartosz Perepeczo Oct 31 '17 at 17:25
  • You have my approach in my answer. I hope that it helps. – xunilk Oct 31 '17 at 23:37
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.

enter image description here

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.

enter image description here

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).

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.