0

I have a TIFF file from MODIS MCD64A1 (burned area) which has a quality flag band with documentation like so:

enter image description here

I want to just select where bit 0 and bit 1 = 1 or where the flags are a land grid cell and also where there was sufficient data. To do this I am first reading in my TIFF file:

library(terra)

#read in quality flag tif
qflags = rast('mcd64a1_flag.tif')

#get unique values in qflags
all_uniques = unique(values(qflags, mat = F))

Now I am looping through the unique decimal values and converting to binary:

#loop through and convert the values to the 8 digit binary code but ignore NaN and negatives as that is bad data already.
#function to convert to 8 bit binary
toBinary <- function(n){
     paste0(as.integer(rev(intToBits(n)[1:8])), collapse = "")
}

#empty vector to store good values
all_good = c()

#start loop

for(i in all_uniques){
  
  #ignore negative and nan
  if(i >= 0 & is.na(i) == FALSE){
    
    #convert the decimal value to binary
    binary = toBinary(i)
    
    #get the last two digits in the binary value
    last_two = as.numeric(substr(binary, nchar(binary) - 1, nchar(binary)))# Extract the last two digits


    #if the last two digits are 11 save the value
    if(last_two == 11){
      
      all_good[[length(all_good) + 1]] = i
    }
    
  }
}

and now I am simply removing all pixels not equal to these good flags:

qflags[!qflags %in% all_good] = NA

but this clearly isn't correct when I look visually and I am not sure where I am going wrong. Is there maybe an easier way or better way to select the good quality flags I need based on the decimal values?

EDIT:

Since bit flags which are water will always be clear in bit flag 1, I guess I really only need to have the bit flag of 1 equal to 1 I think then.

2
  • Why not just use a bitwise operator?
    – Vince
    Sep 10, 2023 at 21:12
  • I guess I am not very familiar with flags stored as bits and I don't know how or haven't seen how to do that. But I am open to anything that will work out. Sep 10, 2023 at 21:26

1 Answer 1

1

This boils down to "select where bit 0 and bit 1 = 1".

Counting up from 0, this happens every fourth integer (ie its true for 3, 7, 11, 15,...). So you can use a modular arithmetic operator, testing if the value mod 4 equals 3.

Reading the data and plotting the quality band (number 3) (I've got an HDF from a NASA Earthdata archive):

mcd = rast("./MCD64A1.A2023182.h12v13.061.2023252124638.hdf")
plot(mcd[[3]])

enter image description here

Then the data with bits 0 and 1 set, ignoring all values in other bits, is this:

plot(mcd[[3]] %% 4 == 3)

enter image description here

For more general bitwise-twiddling of rasters, use the base bitw* functions and apply them to the raster data. For example, the quality flags in bits 5 to 7 encode values from 0 to 7 but you need to extract those bits and shift them right (or divide by enough powers of 2). First write a function:

bit57 = function(x){bitwShiftR(bitwAnd(x, 32+64+128), 5)}

Here I've added 32, 64 and 128 to get a mask that is binary 11100000, and that is and-ed with the data. That zeroes everything except bits 5 to 7. Then the result is shifted right five steps so the three bits we are interested in are now coding 1, 2, and 4 in binary.

This can even be simplified, since a right-shift of 5 will throw out the bits we don't care about, leaving bits 5 to 7 in the right place:

bit57 = function(x){bitwShiftR(x, 5)}

but I'll leave the other function as an illustration of bitwise-AND.

That done, apply it. Here's the map of values encoded by bits 5 to 7 in the quality flag of my sample data:

special = app(mcd[[3]], bit57)
plot(special)

enter image description here

Note you can't call the bitwise functions on terra rasters:

> bitwAnd(mcd[[3]], 64+16)
Error in bitwAnd(mcd[[3]], 64 + 16) : 'a' and 'b' must have the same type

but using app seems to work fine.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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