4

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?

4

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.

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.