When working with MOD09GA product, is supposed QC bitmask for high aerosol quantity to return 4 values?

Question

it's the first time I'm working with MODIS data and I am a bit confused about how bits work and how masks are created. After a bit of research I've been able to understand, more or less, how bit filtering works -- using the help found here and MODIS documentation files (1 and 2).

I am using MOD09GA product to create cloud and aerosol free mask data. I built my binary filter based on the information found in the next table (file 2, page 24).

When doing bitwise operations in python and applying the bit filter to the QA band (state_1km), I am not fully sure that what it returns is right, since it generates a mask with 4 integer values (0, 1, 2 and 3) when I was expecting 2 (0 and 1).

I will show this with an example. You can download the .hdf file from here, although you'll need you login credentials.

import gdal
import numpy as np

# Opening MODIS hdf file
 modisImage = gdal.Open("data/MOD09GAMOD09GA.A2008001.h10v05.006.2015169041625.hdf") 

# Getting subdatasets
subData = modisImage.GetSubDatasets()
data = {}
for fname, name in subData:
    newname = name.split("] ")[1]
    without_type = newname.split(" (")[0]
    print(without_type)
    data[without_type] = gdal.Open(fname).ReadAsArray()

# =============================================================================
# # num_observations_1km MODIS_Grid_1km_2D
# # state_1km_1 MODIS_Grid_1km_2D    <-- I want this one
# # SensorZenith_1 MODIS_Grid_1km_2D
# # .
# # .
# # .
# # .
# # q_scan_1 MODIS_Grid_500m_2D
# =============================================================================


# now call the desired band
qc = data["state_1km_1 MODIS_Grid_1km_2D"]

Cloud state, cloud shadow and land/water bit masks are returning a mask with the values expected (0 and 1), but when I apply the same approach to the aerosol quantity I can't understand the mask returned.

# Aerosol quantity information is in the Bitshift 6 and Bitmask 0-3 (first and third columns of the table). In this case, I want to return only high quantity aerosol

qcAerosol = (qc >> 6) & 3

# Checking values returned (more than 2 values)
np.unique(qcAerosol) # array([  0,  1, 2, 3], dtype=uint16)

# Taking a closer look to te values returned
bin(1)  # '0b1'
bin(2) # '0b10'
bin(3) # '0b11'

The qcAerosol mask created is the one in the image. As you can see, it has four colors corresponding to each value.

Each value returned in the mask, maps a binary value in the table:

• 1 corresponds to the bit mask for 'low aerosol quantity'

• 2 corresponds to the bit mask for 'average aerosol quantity'

• 3 corresponds to the bit mask for 'high aerosol quantity'

Since I want a 'high aerosol quantity' mask, I need a last step to create one.

qcAerosol[qcAerosol != 3] = 0
qcAerosol[qcAerosol == 3] = 1

Getting the next mask.

Python qcAerosol

However, as I noted above, bitmasks for cloud state, cloud shadow and land/water only return two values, 0 and 1. I think there could be a problem with my approach, therefore my question is raised: When working with MOD09GA product, is supposed QC bitmask for high aerosol quantity to return 4 values?

Answer 1

However, as I noted above, bitmasks for cloud state, cloud shadow and land/water only return two values, 0 and 1. I think there could be a problem with my approach, therefore my question is raised: When working with MOD09GA product, is supposed QC bitmask for high aerosol quantity to return 4 values?

According to the description of the quality assessment flags, there can be 4 combinations: '00', '01', '10' and '11'. Thus, it should and it does return 4 values as long as they exist in the QA layer.

Details

The quality assessment (binary) flags for MOD09GA Science Data Sets are indeed detailed in MODIS Surface Reflectance User’s Guide (Collection 6). Specifically, the Data product state QA flags are described in Table 13 (the one you copy-pasted in your question).

Breaking down the instruction

qcAerosol = (qc >> 6) & 3

after the = sign in two parts:

1. qc >> 6 which performs a bitwise right shift of all pixel values in qc by 6 positions. To exemplify, let us use the unsigned 16-bit integer 174272. Its binary form is (all code below in Python):

   bin(174272)[2:]
   '101010100011000000'
   

   Shifting this bit-pattern to the right by 6 positions, will result in

   bin(174272 >> 6)[2:]
   '101010100011'
   

   which is the integer

   int(bin(174272 >> 6)[2:], 2)
   2723
   

   which, by the way, is the same as

   174272 / 2**6
   2723
   

   Theory says we multiply if we shift to the left!)

2. & 3 performs a bitwise AND of all shifted qc values with 3 which is the integer form for the binary string '11'. The 'AND' operation will compare the first two bits of each (now shifted) qc pixel values with '11'.

   Continuing the example with the integer 174272, the second step is:

   int(bin(174272 >> 6)[2:], 2) & 3
   3
   

   This is expected as we compare the first two bits (reading from right to left) '101010100011' with '11' which in turn will result to '11' which is the binary form of the integer 3.

   For the sake of completeness, the possible comparisons are:

   • '00' AND '11' results in '00' or 0
   • '01' AND '11' results in '01' or 1
   • '10' AND '11' results in '10' or 2
   • '11' AND '11' results in '11' or 3

   Thus, np.unique(qcAerosol) can return at most 4 different integer values (0 or 1 or 2 or 3). It may, however, return less or only one according to the recorded state of the pixels for bits 7 and 6.

As for the other data quality states,

• "cloud state" is a double bit, thus 2^2 = 4 states can be encoded
• "cloud shadow" is a single bit, thus only 2^1 = 2 states can be encoded (0 or 1)
• "land/water" is a triple bit, thus the possible combinations are 2^3 = 8: '000' or 0, '001' or 1, '010' or 2, '011' or 3, '100' or 4, '101' or 5, '110' or 6, '111' or 7.

Answer 2

For anyone struggling with bitpacked mask values in MODIS (like I was), I made a package to deal with them more easily. Details are here https://pypi.org/project/unpackqa/

Code to pull the aerosol quantity from the above MOD09 file is as follows

from osgeo import gdal
from matplotlib import pyplot as plt
import unpackqa

# Define how qa/qc flags are from Table 13
mod09ga_spec = {'flag_info':{'cloud_state':[0,1],
                             'cloud_shadow':[2],
                             'land_water_flag':[3,4,5],
                             'aerosol_quantity':[6,7],
                             'cirrus_detected':[8,9],
                             'internal_cloud_algorithm_flag':[10],
                             'internal_fire_algorithm_flag':[11],
                             'MOD35_snow_ice_flag':[12],
                             'pixel_adjacent_to_cloud':[13],
                             'salt_pan':[14]},
                 'max_value' : 65535,
                 'num_bits'  : 16}

# Read directly the "state_1km_1" band. This string is derived from GetSubDatasets() in gdal
qc = gdal.Open('HDF4_EOS:EOS_GRID:"MOD09GA.A2008001.h10v05.006.2015169041625.hdf":MODIS_Grid_1km_2D:state_1km_1').ReadAsArray()

# unpack flags into a dictionary.
flags = unpackqa.unpack_to_dict(qc, product=mod09ga_spec)

unique_values = [0,1,2,3]
plt.imshow(flags['aerosol_quantity'], interpolation='none')
plt.colorbar(ticks=unique_values,values=unique_values)
plt.show()

unpackqa.unpack_to_dict() returns a dictionary where keys are the flag names defined in mod09ga_spec, and each key value is an array the same shape as the qc raster.

flags.keys()
dict_keys(['cloud_state', 'cloud_shadow', 'land_water_flag', 'aerosol_quantity', 'cirrus_detected', 'internal_cloud_algorithm_flag', 'internal_fire_algorithm_flag', 'MOD35_snow_ice_flag', 'pixel_adjacent_to_cloud', 'salt_pan'])

flags['cloud_shadow'].shape
(1200,1200)