Mean of NDVI rasters across time point in QGIS after removing bad quality pixels

Question

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 Cell stats accumulate

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

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
done

# 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

HTH