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 vegetation 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 MODIS Vegetation Index Products (NDVI and EVI), which are two primary vegetation layers; the Normalized Difference Vegetation Index and the Enhanced Vegetation Index. This product will contain data produced on 16-day periods as well as monthly temporal averaged data, in 3 spatial resolutions (250m, 500m and 1km).
The product IDs for each of these products can also be found on the data product page.
This tutorial will use the ‘Vegetation Indices 16-Day L3 Global 250’ product with the product IDs MOD13Q1(Terra Product ID), and MYD13Q1(Aqua Product ID), but will be represented by M*D13Q1 - the second character is replaced by an asterix(*) to identify both Terra and Aqua.
The product layers (original MODIS layers, quality layers and spectral indexes) available for a given product can be retrieved using the following function.
## $prodname ##  "Vegetation Indexes_16Days_250m (M*D13Q1)" ## ## $bandnames ##  "NDVI" "EVI" "VI_QA" "b1_Red" "b2_NIR" "b3_Blue" ##  "b7_SWIR" "View_Zen" "Sun_Zen" "Rel_Az" "DOY" "Rely" ## ## $bandfullnames ##  "16 day NDVI average" "16 day EVI average" ##  "VI quality indicators" "Surface Reflectance Band 1" ##  "Surface Reflectance Band 2" "Surface Reflectance Band 3" ##  "Surface Reflectance Band 7" "View zenith angle of VI pixel" ##  "Sun zenith angle of VI pixel" "Relative azimuth angle of VI pixel" ##  "Day of year of VI pixel" "Quality reliability of VI pixel" ## ## $quality_bandnames ##  "QA_qual" "QA_usef" "QA_aer" "QA_adj_cld" "QA_BRDF" ##  "QA_mix_cld" "QA_land_wat" "QA_snow_ice" "QA_shd" ## ## $quality_fullnames ##  "VI Quality" ##  "VI usefulness" ##  "Aerosol quantity" ##  "Adjacent cloud detected" ##  "Atmosphere BRDF correction performed" ##  "Mixed Clouds" ##  "Land/Water Flag" ##  "Possible snow/ice" ##  "Possible shadow" ## ## $indexes_bandnames ##  "SR" "NDFI" "NDII7" "SAVI" ## ## $indexes_fullnames ##  "Simple Ratio (NIR/Red)" ##  "Flood Index (Red-SWIR2)/(Red+SWIR2)" ##  "NDII7 (NIR-SWIR2)/(NIR+SWIR2)" ##  "SAVI (NIR-Red)/(NIR+Red+0.5)*(1+0.5)"
Note how the
\$bandfullnames define each of the
\$quality_fullnames define the
\$quality_bandnames, and the
\$indexes_fullnames define 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 NDVI (Normalized Difference Vegetation Index) in Mongolia, first we download the boundary of Mongolia 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 Mongolia map_boundary <- geoboundaries("Mongolia") # Defining filepath to save downloaded spatial file spatial_filepath <- "VegetationData/mongolia.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 NDVI data.
To download data in Mongolia, we use the boundary of Mongolia 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 vegetation data can be categorised as spatio-temporal data, the
end_date arguments define the period for which we want the data to be downloaded. Here we use the same date for both the start and end date since we want to download the data for a single date.
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 = "VegetationData", out_folder_mod = "VegetationData", selprod = "Vegetation_Indexes_16Days_1Km (M*D13A2)", bandsel = "NDVI", user = "mstp_test" , password = "MSTP_test_01", start_date = "2020.06.01", end_date = "2020.06.01", 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 NDVI (Normalized Difference Vegetation Index) of Mongolia.
# install.packages(c("rgeoboundaries", "sf, "raster", "here", "ggplot2", "viridis")) library(rgeoboundaries) library(sf) library(raster) library(here) library(ggplot2) library(viridis) # Downloading the boundary of Mongolia map_boundary <- geoboundaries("Mongolia") # Reading in the downloaded NDVI raster data NDVI_raster <- raster(here::here("VegetationData/mongolia/VI_16Days_1Km_v6/NDVI/MYD13A2_NDVI_2020_153.tif")) # Transforming the data NDVI_raster <- projectRaster(NDVI_raster, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") # Cropping the data NDVI_raster <- raster::mask(NDVI_raster, as_Spatial(map_boundary)) # Dividing values by 10000 to have NDVI values between -1 and 1 gain(NDVI_raster)=0.0001 # Converting the raster object into a dataframe NDVI_df <- as.data.frame(NDVI_raster, xy = TRUE, na.rm = TRUE) rownames(NDVI_df) <- c() # Visualising using ggplot2 ggplot() + geom_raster(data = NDVI_df, aes(x = x, y = y, fill = MYD13A2_NDVI_2020_153)) + geom_sf(data = map_boundary, inherit.aes = FALSE, fill = NA) + scale_fill_viridis(name = "NDVI") + labs(title = "NDVI (Normalized Difference Vegetation Index) in Mongolia", subtitle = "01-06-2020", x = "Longitude", y = "Latitude")
The following example uses the
leaflet() function from the
leaflet package to visualise the downloaded NDVI (Normalized Difference Vegetation Index) of Mongolia.
#install.packages("leaflet") library(leaflet) # Defining color palette pal <- colorNumeric(c("#440154FF", "#238A8DFF", "#FDE725FF"), values(NDVI_raster), na.color = "transparent") # Visualising using leaflet leaflet() %>% addTiles() %>% addRasterImage(NDVI_raster, colors = pal) %>% addLegend(pal = pal, values = values(NDVI_raster), title = "NDVI")