3

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.

6

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:

enter image description here

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 i10 format, where bit flags are in i2. 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)

plot

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)

enter image description here

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)

enter image description here

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)

enter image description here


Reference RGB image:

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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