MOD09A1 QC layer

Question

Does anyone know how to handle MOD09A1 Quality flags (QC) in in R? I was trying to use raster package but unable to unlock the QC band.

Answer 1

This is a good question, naturally isn't obvious how to handle this kind of data. For example, the product State QA has 16 bits and there are single, double and triple bit flags combinations:

For example, a bit flag combination as 0-0-0-0-0-0-00-01-001-0-00 means:

• 00 cloud state: clear
• 0 cloud shadow: no
• 001 land/water flag: land
• 01 aerosol quality: low
• 00 cirrus detected: no
• 0 internal cloud algorithm flag: no
• 0 internal fire algorithm flag: no
• 0 MOD35 snow/ice flag: no
• 0 Pixel is adjacent to cloud: no
• 0 BRDF correction performed: no
• 0 internal snow algorithm flag: no

Decoding bit flags like this is really easy, but QA band came in a format, where bit flags are in . So, we need to transform a base 10 number to a base 2 number.

---

R procedure

I'll use a real example for this purpose:

library(raster)

qflags <- raster('~/Downloads/2269893769/MOD09A1_A2017209_h11v11_006_2017218044841_MOD_Grid_500m_Surface_Reflectance_sur_refl_state_500m_deb6f7af.tif')

plot(qflags)

You can convert integers numbers to binary with intToBits() function, but I'll do it from zero creating a function that converts to binary with other options for this purpose:

BitFlagDecoder <- function(x,bits=16,obj=15){
  i <- 1
  result <- rep(0,bits)
  if(length(obj)>1){
    ind <- max(obj):min(obj)+1 # bit index
  }else{
     ind <- obj + 1 
    }
  if(is.na(x)){
    result[1:bits] <- NA
  }else{
    while(x > 0){
      result[i]  <- x %% 2
      x <- x %/% 2
      i <- i + 1
    }
  }
  return(as.numeric(paste(result[ind],collapse = '')))
}

Where bits is the number of bits of the scene and obj is the objective zone or Bit No. (see table above). Just calculate new values and add them to a new raster:

newValues <- sapply(values(qflags),FUN= BitFlagDecoder)

r <- setValues(qflags,newValues)

plot(r)

Also, you can create dictionaries to plot the raster as categorical:

dictionary <-  data.frame(Bit=c(0,1,10,11,100,101,110,111),
                          Class=c('shallow ocean','land',
                                  'ocean coastlines and lake shorelines',
                                  'shallow inland water',
                                  'ephemeral water','deep inland water',
                                  'continental/moderate ocen',
                                  'deep ocean'))

rClass <- as.factor(r)
lvls <- levels(rClass)[[1]]
lvls[['BitFlags']] <- dictionary[match(lvls[,1],dictionary$Bit),'Class']
levels(rClass) <- lvls

library(rasterVis)

levelplot(rClass)

And finally, testing for snow (Bit No. 15, the last bit):

BitFlagDecoder <- function(x,bits=16,obj=15){
  i <- 1
  result <- rep(0,bits)
  if(length(obj)>1){
    ind <- max(obj):min(obj)+1 # bit index
  }else{
     ind <- obj + 1 
    }
  if(is.na(x)){
    result[1:bits] <- NA
  }else{
    while(x > 0){
      result[i]  <- x %% 2
      x <- x %/% 2
      i <- i + 1
    }
  }
  return(as.numeric(paste(result[ind],collapse = '')))
}

newValues <- sapply(values(qflags),FUN= BitFlagDecoder)

r <- setValues(qflags,newValues)

rClass <- as.factor(r)
lvls <- levels(rClass)[[1]]
lvls[['BitFlags']] <- c('No snow','Snow')
levels(rClass) <- lvls

levelplot(rClass)

---

Reference RGB image: