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

Question

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?!

Answer 1

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

Answer 2

I am new to this too so I will not be able to help with prepackages solutions but from a programming perspective we have,
6000 binary rasters as independent files (species exists in pixel or not)
(what dimensions X & Y?)

each file contains an identical but irregular area of interest
( what size bounding box for area of interest?)
( are number of pixels inside bounding box but outside area of interest smallish?)

crop the bounding box out of all rasters
reform each bounding box as a bit vector (use presence=1 absence=0)
stack your 6k bit vectors together and sum in any reasonable language
reform the resulting vector of totals back into a rectangle region. insert the new raster back into its context.

you could get fancy and make a mask (another bit vector) of pixels not to sum over but I would do the easiest first.

I don't know what format your rasters are in now
but converted to netcdf there are tools like nco (netcdf operators)
and ncl (ncar command language) which have many tools to efficiently
manipulate and process rasters.

cheers