I have several NDVI and EVI rasters from different time points across a study period. These are of a region in Burkina Faso, and during the rainy season pixel reliability can be low. I have first masked each raster using the pixel reliability raster so that each NDVI and EVI raster has only the pixels that had a quality score of 0 or 1 (scores are 0, 1 or 3 with 3 being less reliable). I used GRASS r.mask.rast to do this. To get NDVI/EVI across the region I need to combine 6 layers for each index and get the mean values for each pixel.

I have tried to do this using the QGIS raster calculator (my input screenshot 1), cell statistics (my input screenshot 2) and GRASS r.series.accumulate (my input screenshot 3) but I get quite different results for each. The calculator gives me almost no actual pixels back. The cell statistics gives me all the pixels but they look a lot paler (lower values) for the mean than I expected, as does the max value. This might be accurate and I might be being stupid. For the accumulate I get something more like what i would expect with all the pixels, but I do get several warning messages when I run it (the log closes before I can copy them to paste!). All of my rasters are set to ESPG 32630 UTM zone 30N

calculator 1 Cell stats 2 accumulate accumulate

If none of these things is the correct process, is there a better process?

1 Answer 1


Since you started with GRASS GIS, why not do the whole process in using GRASS? There is a module r.series ( manual ) that can calculate any one of several statistics on a pixel-by-pixel basis from a series of rasters. This module also ignores NULL values (which is probably your problem above).

If you have six EVI maps and six quality masks as follows: evi_1, evi_2,...evi_6 and evi_qa_1, evi_qa_2, ... evi_qa_6 then you could try the following loop in the GRASS command shell:

# *Always* set computational region at start
g.region -ap rast=evi_1  # assumes all maps are overlapping
# First loop thru all vegetation index maps and prepare the masked maps
for i in 1 2 3 4 5 6; do
  # set the mask
  r.mask evi_qa_${i} maskcats=1,2  # only use the values 1,2 for the mask
  # create a new, masked veg index map
  r.mapcalc "evi_masked_${i} = evi_${i}"
  # remove the mask for the next loop 
  r.mask -r

# Now build the list of masked veg index maps and then create an averaged series
MAPLIST=`g.list type=raster pattern="evi_masked*" separator=comma`
r.series input=${MAPLIST} output=evi_average method=average


  • Hi Micha, thank you very much for this, which I will try. I have one query - I notice that you first set a mask, and then use the calculating function, before removing the mask. How does this compare to when using the tools separately? I.e., I will have to give a step-by-step of this part during one of my group meetings, and it might be handy to know how to do it both through the shell and through the point-and-click avenues. Again, thank you very much!
    – roztaylor
    Jun 27, 2022 at 18:26
  • Not sure I understand your question. When you set a mask in GRASS, all subsequent actions use that mask. So every result for whatever module you use will have values only within the mask. Since your quality rasters for the VI data will probably be different each time point, then you need to set the mask, run r.mapcalc to get the masked EVI values, then remove the mask to prepare for the next time point.
    – Micha
    Jun 28, 2022 at 20:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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