Setting pixels of a all raster layers in a raster stack to NA based on the pixel values of corresponding raster layers of another raster stack

Question

I have MODIS NDVI product (MOD13A2) for 10 years in HDF file format. Here i have 230 of such HDFs. For a single HDF, I want to set NDVI pixel values to NA with respect to Quality values in another Quality band. I got 230 NDVI HDFs in one folder and corresponding 230 Quality HDFs in another folder. E.g. for quality band pixel values of 1,2 and 3 i want to set NDVI values to NA.

How can i perform this work on all 230 hdfs.

I am using R and here is the code i tried:

library(raster)
# read all files with .hdf extension in R    
NDVI_data<-list.files(path, pattern= "\\.hdf$", all.files=FALSE, full.names=TRUE,recursive=FALSE) 
NDVI_reliability_data<-list.files(path, pattern= "\\.hdf$", all.files=FALSE, full.names=TRUE,recursive=FALSE) 
s_NDVI <- stack(NDVI_data)
s_NDVI_reliability = stack(NDVI_reliability_data)

I want one on one comparison of 1,2,...230 layers of s_NDVI to the corresponding 1,2,...230 of s_NDVI_reliability so that pixels of s_NDVI layers set to NA based on pixels values of s_NDVI_reliability layers.

Answer 1

If you have two rasters of the same shape you can replace conditionally much like you would a vector in R. For example:

> ndvi = raster(matrix(1:12,3,4))
> qual = raster(matrix(sample(1:12),3,4))

To set ndvi to NA where qual <= 3, do:

> ndvi[qual <= 3] = NA

and to check, here's the corresponding values of ndvi and qual after the replacement - all the NAs are in the right place:

> cbind(as.data.frame(ndvi), as.data.frame(qual))
   layer layer
1      1     7
2      4     8
3     NA     3
4     10    12
5      2     4
6     NA     2
7      8    11
8     NA     1
9      3     5
10     6     9
11     9     6
12    12    10
>