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?
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.
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.
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")
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.