[GIS] How to sum values of overlapped pixels in thousands of rasters layers and clip it using a polygon in QGIS

qgisrastersaga

I have a lot of rasters files (more than 6000), each one represents the presence and absence of a species (1 for presence and NoData for absence). The set of rasters covers the whole Latin America and they can overlap eventually. I also have a vector polygon layer which bounds an area in which I am interested in. As a final goal, I would like to sum all the species in each grid in my area of study.
Any ideas on how to achieve this?
I have tried to sum all raster using Grid Sum in Saga and in order to clip it afterwards. Problem is that the output has only NaN and zeros values – which does not make sense… I used a sample of 100 rasters, in order to test it. And it lasts more than two hours! Any ideas on how to do it faster?!

Best Answer

Following on Vince's comment: you can remap your species to IDs that are powers of 2, as below, represented as integers and binary.

Species 1 --> 1 --> 00000001
Species 2 --> 2 --> 00000010
Species 3 --> 4 --> 00000100

You can see that each species is identified by a 1 in a unique slot in the bit sequence. So for each species, you'd have a raster with cells of either zero (no presence) or, say, 4 (species 3 is present). In a 32-bit integer, you have 32 of these slots.

When you add these layers together, you get a sum for each pixel that indicates which species are present.

Raster 1 --> 000000001
Raster 2 --> 000000010
Raster 3 --> 000000100
             ---------
Output   --> 000000111 --> 6

In this instance, the value 6 tells you that species 1, 2 and 3 are present.

You can also use other operators such as bitwise AND not NOT to see cells where two or more specific species are present, or cells where certain species are not present together.

If you process your rasters this way, you can add them incrementally using a scriptable raster calculator like gdal_calc.py, adding each raster to the the sum of all previous rasters (though you'd probably have to turn all your rasters into VRTs first so that they have equal extents, etc.)

In QGIS you'll just add them all to the summation in the raster calculator.

Edit, based on the answer to Jens' comment and more info:

I'm answering the wrong question above. If you want only to sum the number of species in a cell and each cell-with-a-species is 1, you just add all the rasters together. To convert the NaN cells to zero, you can use the numpy function nan_to_num, which converts NaNs to zeroes. So you add like this:

nan_to_num("raster1@1") + nan_to_num("raster2@1")
Related Question