Quality control

Not all pixels in a MODIS image are not suitable for analysis of land areas. For example, pixel quality can be negatively affected by clouds or other atmospheric conditions. For that reason, each image needs to be pre-processed to remove the values that cannot be used.

Each MODIS file contains quality assurance (QA) information that can be used to identify the pixels to remove. The QA information is stored in a somewhat complicated bit-encoding format. This allows for very efficient storage, but it makes it much harder to use.

A single bit can represent two values — 0 (no) or 1 (yes); a combination of two bits can represent 2² = 4 values (00, 01, 10, and 11), and with three bits you can store 2³ = 8 values. This video by USGS explains the QA encoding.

The table below shows State QA description (16-bit) for 500 m, 1 km and lower resolution MODIS surface reflectance products. For product specific QA attribute, follow the MODIS Surface Reflectance User’s Guide.

bit	variable	value	description
1-2	cloud state	00	clear
		01	cloudy
		10	mixed
		11	not set, assumed clear
3	cloud shadow	1	yes
		0	no
4-6	land/water flag	000	shallow ocean
		001	land
		010	ocean coastlines and lake shorelines
		011	shallow inland water
		100	ephemeral water
		101	deep inland water
		110	continental/moderate ocean
		111	deep ocean
7-8	aerosol quantity	00	climatology
		01	low
		10	average
		11	high
9-10	cirrus detected	00	none
		01	small
		10	average
		11	high
11	cloud flag	1	cloud
		0	no cloud
12	fire flag	1	fire
		0	no fire
13	snow/ice flag	1	yes
		0	no
14	is adjacent to cloud	1	yes
		0	no
15	Salt pan	1	yes
		0	no
16	Snow mask	1	yes
		0	no

Note that in R, the first bit will be referred as “1”, whereas in many other languages (e.g. Python) it will follow the values as shown in table.

To interpret the pixel level QA values, we need to convert them from decimal to binary format. The luna package offers a function for this conversion and create a mask from the QA band. The user need to specify a matrix (“qabits”) with the start and end of the quality assessment (QA) bits considered, and specify a list (“reject”) with the values to be rejected (in the image) matching the rows in qabits. Following the table above, we will define the “reject” values to exclude pixels affected by cloud and cloud shadow.

from <- c(1,3,11,14)
to   <- c(2,3,11,14)
reject <- c("01,10", "1", "1", "1")
qa_bits <- cbind(from, to, reject)
qa_bits
##      from to   reject
## [1,] "1"  "2"  "01,10"
## [2,] "3"  "3"  "1"
## [3,] "11" "11" "1"
## [4,] "14" "14" "1"

Pixels with bits 1 and 2 with values “01” or “10” will be rejected. All other combinations (“00” and “11” in this case) are not rejected.

We use the downloaded MODIS file that, in a previouse step, we saved in the datadir directory.

library(terra)
datadir <- file.path(dirname(tempdir()), "_modis")
mf <- file.path(datadir, "MOD09A1.A2009361.h21v08.061.2021149144347.hdf")
r <- rast(mf)

Generate the quality mask. We will use band 12 sur_refl_state_500m that has the quality data.

qc <- r[[12]]
plot(qc, main = "Quality")

The luna package has a modis_mask method to create a mask from the quality band and the parameters defined above.

library(luna)
##
## Attaching package: 'luna'
## The following object is masked from 'package:base':
##
##     version
quality_mask <- modis_mask(qc, 16, qa_bits)
plot(quality_mask, main="Quality mask")

The plot shows the pixels we want to retain. Now that we have the quality mask, we can apply it to all the bands. It is always useful to visually check the result of the masking as many times pixels could be wrongly flagged as “poor” quality. Residual noises can be filtered in subsequent gap-filling and smoothing operations.

rmask <- mask(r, quality_mask)

And we can plot the results, here as a “false color composite” (NIR:Red:Green)

plotRGB(rmask, r = 2, g = 1, b = 4, main='False color composite', stretch="lin")

Finally we save the result after cloud masking.