Accessing the Quality Band for MCD12Q1 with MODISTools GetBands

Question

I am using MODISTools (https://github.com/seantuck12/MODISTools) to process MCD12Q1 data, which, according to the LP DAAC website (https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mcd12q1), has 16 bands. I believe one of them is a quality band (Land Cover QC). However, when I use the MODISSubsets function to download MCD12Q1 and the GetBands function to view the bands, only 15 bands are returned (without the Land Cover QC band). Without this band, I can't use the MODISSummaries function. Does anyone know how to work around this issue?

Answer 1

It is easy to extract needed subdatasets from HDF with GDAL command line tools. I guess you might want to do everything in R so here is how they can be run from it. In R there is a package called gdalUtils it provides wrappers for GDAL tools. You'll need gdal_translate for data conversion.

Example:

library(gdalUtils)
hdf4_dataset <- file.choose() #select your HDF
gdal_translate(hdf4_dataset,"modis_qa.tif",sd_index=1)

You can then read it with another R package rgdal:

library(rgdal)
d <- file.choose() #select your saved GeoTIFF
x = new("GDALReadOnlyDataset", test)

Now manipulate x to your liking it is just a data.frame.

You might wonder why not use rgdal only. The problem is that it doesn't read HDF format which is used to distribute MODIS data.

Answer 2

This can easily be achieved using the MODIS package which is currently hosted on GitHub. Install the package via

library(devtools)
install_github("MatMatt/MODIS")
library(MODIS)

and make sure to adjust the local output filepaths via MODISoptions. A more detailed description of the installation procedure can e.g. be found here.

Next, download a desired MCD12Q1 tile using

getHdf("MCD12Q1", begin = "2013-01-01", 
       tileH = 21, tileV = 9)

and, finally, list the scientific datasets (SDS) included therein using

fls <- list.files(getOption("MODIS_localArcPath"), 
                  pattern = "MCD12Q1.A2013001.*.hdf",
                  recursive = TRUE, full.names = TRUE)

getSds(fls)

As seen from the console output, the 11th SDS is the desired 'Land Cover Type QC' layer, which can subsequently be extracted using

runGdal("MCD12Q1", begin = "2013-01-01", 
        tileH = 21, tileV = 9, 
        SDSstring = "0000000000100000")

Note that runGdal is also a wrapper function around getHdf, meaning that if the desired MCD12Q1 tile has not been previously downloaded, runGdal will automatically perform the download and extract the desired SDS layers.