Using Quality Control process for MODIS Leaf Area Index?

Question

Anyone can help in explaining how I can apply Quality Control (QC) for MODIS Leaf Area Index (LAI) to extract the valid pixels of LAI using ArcMap?

Answer 1

I assume that you are referring to one of the following MODIS products: MYD15A2, MOD15A2, MCD15A2 or MCD15A3.

The quality control layer for these files is based on using each of the bits in the value as separate flags. ArcMap doesn't know this, and interprets it as a standard 8-bit image (values from 0 to 255). The intended approach would be to use LDOPE to mask the data, but it can be a somewhat daunting task to get a handle on.

Another option is to take a very direct approach, and use reclassify in ArcMap to establish a simplified QC layer. Looking at pages 3 and 4 in this PDF will give you a list of QC flags in the product. These can be translated using this converter.

It should be noted that the order of bits in the MODIS QC layer is from right to left. As such, the simplest usable QC mask you can make from the MODIS QC product would be to only use data where the MODIS QC layer is 0 or 2 (meaning Good Quality & either from the Aqua or Terra satellite).

Answer 2

This is an R solution (not ArcGIS) but for reference I though that I would go ahead and post this code soultion. You can parse the value(s) for a given bit(s) based on the 8-bit integer values of the QC raster. Quality flags are stored in each bit, or combination of bits and stored as an integer. For example the bit values for 100 are "0 1 1 0 0 1 0 0" (reversed for right to left storage). If you are interested in bit five then you would just grab the value in the fifth position of the bit string "0 1 1 0 0 1 0 0" which would be zero, indicating the value of the flag.

This is a function to parse the bits based on an integer. This is in the development version of the spatialEco package so, if you use this function please reference the package. I would highly recommend not changing the default values for depth and order. For the MODIS QC data the bits are reversed (right to left) so, the default will parse the bit values correctly. The depth of the bits are also 8, however there could be other applications where this value may need to be larger (eg., 32 for 16 bit integer data).

parse.bits <- function(x, bit, depth=8, order = c("reverse", "none") ) {
  position <- depth - bit
  if(order[1] == "reverse") sort.order = rev(1:depth) else sort.order = 1:depth  
    b <- as.integer(intToBits(x)[sort.order])[depth:1][position]
    if(length(bit) > 1) b <- paste0(b, collapse=" ")
  return(b)
}

Here is an applied example using QC values from the Harmonized Landsat Sentinel-2 product (handbook with QC values, pp 21). First we create some dummy data representing a data raster (r) and a QC band (qc).

library(raster)
r <- raster(nrow=100, ncol=100)
  r[] <- round(runif(ncell(r), 0,1)) 
qc <- raster(nrow=100, ncol=100)
  qc[] <- round(runif(ncell(qc), 64,234)) 

Now, we create a data.frame with bit values from QC table. Please note that you would have to modify this data.frame based on the QC values for any given product. These are documented in the users handbook(s). I would point out that most of these flags are only calling one bit however, if you look at the aerosol flags you will see that they are using bits 7 and 6, yielding 4 possible values, and yes order matters. Some bits are flipped so, pay special attention to how they are defined in the users handbook. Many of the MODIS products QC values utilize multiple bits thus returning more than two flag values.

( qc_bits <- data.frame(int=0:255, 
                      cloud = unlist(lapply(0:255, FUN=parse.bits, bit=1)),
                      shadow = unlist(lapply(0:255, FUN=parse.bits, bit=3)),
                      acloud = unlist(lapply(0:255, FUN=parse.bits, bit=2)),
                      cirrus = unlist(lapply(0:255, FUN=parse.bits, bit=0)),
                      aerosol = unlist(lapply(0:255, FUN=parse.bits, bit=c(7,6)))) )

Now we can query the results to create a vector of integer values indicating what to mask cloud is in bit 1 and shadow bit 3. The possible values in these two bits are 0 (no) and 1 (yes). The result is an integer vector that can be used to set NA values in the QC band and then used to mask the data.

m <- sort(unique(qc_bits[c(which(qc_bits$cloud == 1),
                           which(qc_bits$shadow == 1)
                           ),]$int))

Finally we can use the queried integer values to mask image with QA band

qc[qc %in% m] <- NA
r <- mask(r, qc)