MODIS (Moderate Resolution Imaging Spectroradiometer)

MODIS is an instrument aboard the Terra and Aqua satellites, which orbits the entire Earth every 1-2 days, acquiring data at different spatial resolutions. The data acquired by MODIS describes features of the land, oceans and the atmosphere. A complete list of MODIS data products can be found on the MODIS website.

Downloading Land Cover data using MODIStsp

MODIStsp is an R package for downloading and preprocessing time series of raster data from MODIS data products. The package’s name is an acronym for ‘MODIS Time Series Processing’. This tutorial will focus on downloading and visualising land cover data, but the same process can be followed with other MODIS data products as well.

Installing MODIStsp

The MODIStsp package can be downloaded from CRAN as follows.


Identifying the MODIS data product

The first step of downloading data is to identify which MODIS data product to use.

This tutorial will use the Land Cover Products, which uses different classification schemes to identify types of land cover. The primary land cover scheme identifies 17 classes defined by the International Geosphere-Biosphere Programme (IGBP) as detailed in their user guide.This product contains annually averaged data, collected bi-annually allowing for hemispheric differences in the growing seasons, and enabling the product to capture two growth cycles if necessary.

The product IDs for each of these products can also be found on the data product page.

This tutorial will use the ‘Land Cover Type Yearly L3 Global 500m’ product with the product IDs MCD12Q1 (Combined Aqua and Terra Product ID)

Retrieving MODIS layers for a product

The product layers (original MODIS layers, quality layers and spectral indexes) available for a given product can be retrieved using the following function.

## $prodname
## [1] "LandCover_Type_Yearly_500m (MCD12Q1)"
## $bandnames
##  [1] "LC1"                 "LC2"                 "LC3"                
##  [4] "LC4"                 "LC5"                 "LC_Prop1"           
##  [7] "LC_Prop2"            "LC_Prop3"            "LC_Prop1_Assessment"
## [10] "LC_Prop2_Assessment" "LC_Prop3_Assessment" "LC_QC"              
## [13] "LC_LW"              
## $bandfullnames
##  [1] "Land Cover Type 1 (IGBP)*"                                      
##  [2] "Land Cover Type 2 (UMD)*"                                       
##  [3] "Land Cover Type 3 (LAI/fPAR)*"                                  
##  [4] "Land Cover Type 4 (NPP/BGC)*"                                   
##  [5] "Land Cover Type 5: Annual Plant Functional Types classification"
##  [6] "FAO-Land Cover Classification System 1 (LCCS1) land cover layer"
##  [7] "FAO-LCCS2 land use layer"                                       
##  [8] "FAO-LCCS3 surface hydrology layer"                              
##  [9] "LCCS1 land cover layer confidence"                              
## [10] "LCCS2 land use layer confidence"                                
## [11] "LCCS3 surface hydrology layer confidence"                       
## [12] "Land Cover QC"                                                  
## [13] "Binary land/water mask derived from MOD44W"                     
## $quality_bandnames
## $quality_fullnames
## $indexes_bandnames
## $indexes_fullnames

Note how the \$bandfullnames define each of the \$bandnames.

The MODIStsp() function

MODIStsp() is the main function of the MODIStsp package, and allows us to download MODIS data products. While this is a very comprehensive function and we only a very few of its arguments in this tutorial, the entire list of arguments can be found in the MODIStsp documentation.

The MODIStsp() function provides two ways of downloading data, namely, through a GUI (interactive) or through an R script (non-interactive). This tutorial will focus on the non-interactive execution.

Downloading IGBP land cover classification data

To download the IGBP land cover classification data in Zimbabwe, first we download the boundary of Zimbabwe with the geoboundaries() function from the rgeoboundaries package and save it on our computer.

# install.packages(c("rgeoboundaries", "sf"))
# Downloading the country boundary of Zimbabwe
map_boundary <- geoboundaries("Zimbabwe")
# Defining filepath to save downloaded spatial file
spatial_filepath <- "LandCoverData/zimbabwe.shp"
# Saving downloaded spatial file on to our computer
st_write(map_boundary, paste0(spatial_filepath))

Then we use the MODIStsp() function to download the IGBP land cover classification data.

To download data in Zimbabwe, we use the boundary of Zimbabwe we downloaded. So in the MODIStsp() function we set the spatmeth argument as “file” and set the spafile argument as the path of the map we saved.

Since IGBP land cover classification data is provided annually and dated on the 1st of January each year, the start-date and end_date arguments should contain the 1st of January of each year for which we want the data to be downloaded.

The start-date and end_date arguments define the period for which we want the data to be downloaded.

Note that this tutorial uses a test username and password. The user and password arguments should be the username and password corresponding to your earthdata credentials.

MODIStsp(gui             = FALSE,
         out_folder      = "LandCoverData",
         out_folder_mod  = "LandCoverData",
         selprod         = "LandCover_Type_Yearly_500m (MCD12Q1)",
         bandsel         = "LC1", 
         user            = "mstp_test" ,
         password        = "MSTP_test_01",
         start_date      = "2019.01.01", 
         end_date        = "2019.12.31", 
         verbose         = FALSE,
         spatmeth        = "file",
         spafile         = spatial_filepath,
         out_format      = "GTiff")

Understanding the downloaded data

The downloaded files are saved in subfolders within the defined output folder.

A separate subfolder is created for each processed original MODIS layer, Quality Indicator or Spectral Index with an image for each processed date. The images will be placed in the following folder structure and named using the following naming convention.


Visualising land cover in Zimbabwe

The following example uses the geom_raster() function from the ggplot2 package to visualise the downloaded IGBP land cover classification in Zimbabwe

# install.packages(c("rgeoboundaries", "sf", "raster", "here", "ggplot2", "viridis"))

# Downloading the boundary of Zimbabwe
map_boundary <- geoboundaries("Zimbabwe")

# Reading in the downloaded landcover raster data
IGBP_raster <- raster(here::here("LandCoverData/zimbabwe/LandCover_Type_Yearly_500m_v6/LC1/MCD12Q1_LC1_2019_001.tif"))

# Transforming data
IGBP_raster <- projectRaster(IGBP_raster, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

# Cropping data
IGBP_raster <- raster::mask(IGBP_raster, as_Spatial(map_boundary))

# Converting the raster object into a dataframe and converting the IGBP classification into a factor
IGBP_df <-, xy = TRUE, na.rm = TRUE) %>%
  mutate(MCD12Q1_LC1_2019_001 = as.factor(round(MCD12Q1_LC1_2019_001)))
rownames(IGBP_df) <- c()
# Renaming IGBP classification levels
levels(IGBP_df$MCD12Q1_LC1_2019_001) <- c( "Evergreen needleleaf forests",
                                           "Evergreen broadleaf forests",
                                           "Deciduous needleleaf forests",
                                           "Deciduous broadleaf forests",
                                           "Mixed forests",
                                           "Closed shrublands",
                                           "Open shrublands",
                                           "Woody savannas",
                                           "Permanent wetlands",
                                           "Urban and built-up lands",
                                           "Cropland/natural vegetation mosaics",
                                           "Snow and ice",
                                           "Water bodies")
# Visualising using ggplot2
ggplot() + 
  geom_raster(data = IGBP_df,
              aes(x = x, y = y, fill = MCD12Q1_LC1_2019_001)) +
  geom_sf(data = map_boundary, inherit.aes = FALSE, fill = NA) +
  scale_fill_viridis(name = "Land Cover Type", discrete = TRUE) +
  labs(title = "Land Cover classification in Zimbabwe",
       subtitle = "01-01-2019 - 31-12-2019",
       x = "Longitude",
       y = "Latitude")