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.
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.
MODIStsp package can be downloaded from CRAN as follows.
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)
The product layers (original MODIS layers, quality layers and spectral indexes) available for a given product can be retrieved using the following function.
## $prodname ##  "LandCover_Type_Yearly_500m (MCD12Q1)" ## ## $bandnames ##  "LC1" "LC2" "LC3" ##  "LC4" "LC5" "LC_Prop1" ##  "LC_Prop2" "LC_Prop3" "LC_Prop1_Assessment" ##  "LC_Prop2_Assessment" "LC_Prop3_Assessment" "LC_QC" ##  "LC_LW" ## ## $bandfullnames ##  "Land Cover Type 1 (IGBP)*" ##  "Land Cover Type 2 (UMD)*" ##  "Land Cover Type 3 (LAI/fPAR)*" ##  "Land Cover Type 4 (NPP/BGC)*" ##  "Land Cover Type 5: Annual Plant Functional Types classification" ##  "FAO-Land Cover Classification System 1 (LCCS1) land cover layer" ##  "FAO-LCCS2 land use layer" ##  "FAO-LCCS3 surface hydrology layer" ##  "LCCS1 land cover layer confidence" ##  "LCCS2 land use layer confidence" ##  "LCCS3 surface hydrology layer confidence" ##  "Land Cover QC" ##  "Binary land/water mask derived from MOD44W" ## ## $quality_bandnames ## NULL ## ## $quality_fullnames ## NULL ## ## $indexes_bandnames ## NULL ## ## $indexes_fullnames ## NULL
Note how the
\$bandfullnames define each of the
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.
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.
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")) library(rgeoboundaries) library(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
end_date arguments should contain the 1st of January of each year for which we want the data to be downloaded.
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
password arguments should be the username and password corresponding to your earthdata credentials.
library(MODIStsp) 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")
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.
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")) library(rgeoboundaries) library(sf) library(raster) library(here) library(ggplot2) library(viridis) library(dplyr) # 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 <- as.data.frame(IGBP_raster, 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", "Savannas", "Grasslands", "Permanent wetlands", "Croplands", "Urban and built-up lands", "Cropland/natural vegetation mosaics", "Snow and ice", "Barren", "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")
MODIS website: https://modis.gsfc.nasa.gov/
MODIStsp vignette: https://cran.r-project.org/web/packages/MODIStsp/vignettes/MODIStsp.html
MODIStsp land cover data product page: https://modis.gsfc.nasa.gov/data/dataprod/mod12.php