# MODIS MCD641 quality bits from TIFF file in R

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

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)

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.

• Why not just use a bitwise operator? 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

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]])

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

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

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)

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.