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

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

• How many total species are involved? If it's under 31, you could use a bitwise mask and reclass each data value as a power of 2 (1,2,4,8...) and allow the results to coexist in one raster. Jul 22, 2015 at 15:03
• what is the meaning of summing the rasters? do you want to have the amount of species per pixel? if so, then why do you have -1 for presence? I am sure, recoding all -1 to 1 would help for the Grid Sum in SAGA. clipping is a piece of cake then.
– Jens
Jul 22, 2015 at 19:57
• Yes,Jens, the meaning of summing the raster is to have the amount of species per pixel. And the values for presence is 1 (sorry, my mistake, I edited the question and shift " -" for parenthesis). And Grid Sum is not working, beacause the output is an empty map, with values from NaN and zeros... Jul 24, 2015 at 13:59

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.

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")
``````
• Sorry, Rob, I am new in Qgis and I don't know how to remap rasters... Even if I learn it, is it possible to do that for 6000 rasters??? Jul 24, 2015 at 14:05
• When I say "re-map" I mean in a functional, not cartographic, sense. Like, if you have a species with an integer ID, it would "map onto" an ID that is a power of 2. How many species are we talking about? Are there many rasters for one species, or is each one for a different species? Jul 26, 2015 at 6:07
• Each raster is for one species. I have more than 6000 species (consequently more than 6000 rasters). All of them are the same size, the same number of rows and columns, and the same pixel size. And they overlap. I suspect Grid Sum (SAGA) is not working because of the NoData values... Does it make sense? Maybe if I shift the NoData for zeros? Jul 27, 2015 at 17:22
• I see now. I've answered the wrong question. See the bottom of my answer for more. Jul 27, 2015 at 21:38
• Thank you, very much! Yes, when I changed the NaN to zeroes, Grid Sum worked! Jul 30, 2015 at 13:15

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

• Thank you, very much! I had already solved the problem with Rob suggestion... Jul 30, 2015 at 13:16