In lbusett/MODIStsp: Find, Download and Process MODIS Land Products Data

Introduction

{MODIStsp} is an R package allowing to automatize the creation of time series of rasters derived from MODIS Land Products data. It allows performing several preprocessing steps on MODIS data available within a given time period.

Development of {MODIStsp} started from modifications of the ModisDownload R script by Thomas Hengl [-@Hengl2010], and successive adaptations by Babak Naimi [-@Naimi2014]. The basic functionalities for download and preprocessing of MODIS datasets provided by these scripts were gradually incremented with the aim of:

• developing a stand-alone application allowing to perform several preprocessing steps (e.g., download, mosaicing, reprojection and resize) on all available MODIS land products by exploiting a powerful and user-friendly GUI front-end;
• allowing the creation of time series of both MODIS original layers and additional Quality Indicators (e.g., data acquisition quality, cloud/snow presence, algorithm used for data production, etc. ) extracted from the aggregated bit-field QA layers;
• allowing the automatic calculation and creation of time series of several additional Spectral Indexes starting form MODIS surface reflectance products.

All processing parameters can be easily set with a user-friendly GUI, although non-interactive execution exploiting a previously created Options File is possible. Stand-alone execution outside an R environment is also possible, allowing to use scheduled execution of MODIStsp to automatically update time series related to a MODIS product and extent whenever a new image is available.

Required MODIS HDF files are automatically downloaded from NASA servers and resized, reprojected, resampled and processed according to user's choices. For each desired output layer, outputs are saved as single-band rasters corresponding to each acquisition date available for the selected MODIS product within the specified time period. "R" RasterStack objects with temporal information as well as Virtual raster files (GDAL vrt and ENVI META files) facilitating access to the entire time series can be also created.

Installation

{MODIStsp} requires R v >= 3.6.3.

On Windows

You can install the stable version of {MODIStsp} from CRAN:

install.packages("MODIStsp")

, or the development version (containing the latest improvements and bug fixes) from GitHub:

install.packages("remotes")
library(remotes)
install_github("ropensci/MODIStsp")

On Linux systems

To install {MODIStsp} on Linux, you need to be able to install the {sf} package, which requires several dependencies. See here if you have trouble installing {sf}.

In addition, you need to install dependencies required by the {protolite} package, required by {geojson}. See here for instructions on installing them.

Then, you can install the stable version of {MODIStsp} from CRAN:

install.packages("MODIStsp")

, or the development version (containing the latest improvements and bug fixes) from GitHub:

library(devtools)
install_github("ropensci/MODIStsp")

On Mac OS

To install {MODIStsp} on MacOS, you need to be able to install the {sf} package, which requires several dependencies. See HERE if you have trouble installing {sf}.

Then, you can install the stable version of {MODIStsp} from CRAN:

install.packages("MODIStsp")

, or the development version (containing the latest improvements and bug fixes) from GitHub:

library(devtools)
install_github("ropensci/MODIStsp")

Running the tool in Interactive Mode: the MODIStsp GUI

The easiest way to use {MODIStsp} is to use its powerful GUI (Graphical User Interface) for selection of processing options, and then run the processing.

To open the GUI, load the package and launch the MODIStsp function, with no parameters:

library(MODIStsp)
MODIStsp()

This opens a GUI from which processing options can be specified and eventually saved (or loaded).

The GUI allows selecting all processing options required for the creation of the desired MODIS time series. The main available processing options are described in detail in the following.

  library(png)
  library(grid)
  library(knitr)
  img <- readPNG("GUI_1.PNG")
  grid.raster(img)

See (https://docs.ropensci.org/MODIStsp/articles/interactive_execution.html) for further info and instructions regarding the usage of the GUI.

Non-Interactive Execution from within R

MODIStsp() can be launched in non-interactive mode within an R session by setting the optional GUI parameter to FALSE, and either providing the desired processing argument in the call to the function, or providing a previously saved opts_file specifying the path to a JSON Options file previously saved through the GUI. This allows exploiting {MODIStsp} functionalities within generic R processing scripts.

Specifying the processing parameters in the function call

All processing parameters can be set in the call to MODIStsp(). Mandatory parameters are selprod (specifying the MODIS product), (one of) bandsel, quality_bandsel or indexes_bandsel, (that specify the desired output layers), out_folder and start_date and end_date. user and password are also needed if download_server is not equal to "offline".

The function MODIStsp_get_prodlayers() allows to easily retrieve the names of products and available layers based on product code, such as in:

library(MODIStsp) 
MODIStsp_get_prodlayers("M*D13Q1")

The other parameters are set automatically to default values (see the function reference for details on the different available function arguments).

For example, the following code processes layers NDVI and EVI and quality indicator usefulness of product M*D13Q1, considering both Terra and Aqua platforms, for dates comprised between 2020-06-01 and 2020-06-15 and saves output to R tempdir():

library(MODIStsp) 

# **NOTE** Output files of examples are saved to file.path by setting out_folder to $tempdir.

# --> See name and available layers for product M*D13Q1
# MODIStsp_get_prodlayers("M*D13A2")

# --> Launch the processing
MODIStsp(
  gui             = FALSE, 
  out_folder      = "$tempdir", 
  selprod         = "Vegetation_Indexes_16Days_1Km (M*D13A2)",
  bandsel         = c("EVI", "NDVI"), 
  quality_bandsel = "QA_usef", 
  indexes_bandsel = "SR", 
  user            = "mstp_test" ,
  password        = "MSTP_test_01",
  start_date      = "2020.06.01", 
  end_date        = "2020.06.15", 
  verbose         = FALSE,
  parallel        = FALSE
)

# Outputs are in this case in subfolder "MODIStsp/VI_16Days_1Km_v61" of 
# `base::tempdir()`: 

out_fold <- file.path(tempdir(), "MODIStsp/VI_16Days_1Km_v61/") 
list.files(out_fold)
list.files(file.path(out_fold ,"EVI"))
list.files(file.path(out_fold ,"QA_usef"))

Note that this example, as well as the following ones, is run in single core to follow CRAN policies, by setting parallel = FALSE; users can exploit multicore functionalities skipping to set this argument.

Launching MODIStsp using a saved "Options file"

Alternatively, you can run MODIStsp() without opening the GUI by specifying a previously saved options file:

library(MODIStsp) 

# **NOTE** Output files of examples are saved to file.path(tempdir(), "MODIStsp").

# --> Specify the path to a valid options file saved in advance from MODIStsp GUI 
# Here we use a test json file saved in MODIStsp installation folder which
# downloads and processed 3 MOD13A2 images over the Como Lake (Lombardy, Italy)
# and retrieves NDVI and EVI data, plus the Usefulness Index Quality Indicator.

opts_file <- system.file("testdata/test_MOD13A2.json", package = "MODIStsp")

# --> Launch the processing
MODIStsp(gui = FALSE, opts_file = opts_file, verbose = FALSE, parallel = FALSE)

# Outputs are in this case in subfolder "MODIStsp/VI_16Days_1Km_v61" of 
# tempdir(): 

out_fold <- file.path(tempdir(), "MODIStsp/VI_16Days_1Km_v61") 
list.files(out_fold)
list.files(file.path(out_fold ,"EVI"))

Looping over different Options files

If you need to process different MODIS products, you can prepare beforehand different MODIStsp options files by using the GUI, and then loop over them like this:

opts_files <- c(system.file("testdata/test_MOD13A2.json", package = "MODIStsp"), 
                system.file("testdata/test_MOD10A2.json", package = "MODIStsp"))

for (opts_file in opts_files) {
  MODIStsp(gui = FALSE, opts_file = opts_file, verbose = FALSE, parallel = FALSE)
}

# MOD13A2 ouptuts
out_fold <- file.path(tempdir(), "MODIStsp/VI_16Days_1Km_v61") 
list.files(out_fold)
list.files(file.path(out_fold ,"NDVI"))

# MOD10A2 ouptuts
out_fold <- file.path(tempdir(), "MODIStsp/Surf_Temp_8Days_1Km_v61") 
list.files(out_fold)
list.files(file.path(out_fold ,"LST_Night_1km"))

Specifying the processing parameters using a previously saved options file and overwriting some parameters

Finally, it is possible to both specify a previously saved options file AND setting some parameters in the call to the function. This allows easily performing similar processings, by only updating the required arguments, as in the examples below.

Looping over different spatial extents

Specifying the spafile parameter while setting the spatmeth parameter to "file" overrides for example the output extent of the selected Options File. This allows performing the same preprocessing on different extents using a single Options File. For example:

# Run the tool using the settings previously saved in a specific option file
# and specifying the extent from a spatial file allows to re-use the same
# processing settings to perform download and reprocessing on a different area
library(MODIStsp) 
opts_file    <- system.file("testdata/test_MOD13A2.json", package = "MODIStsp")
spatial_file <- system.file("testdata/lakeshapes/garda_lake.shp", package = "MODIStsp")
MODIStsp(
  gui = FALSE, 
  opts_file = opts_file,
  spatmeth = "file", 
  spafile = spatial_file, 
  verbose = FALSE, 
  parallel = FALSE
)

# --> Create a character array containing a list of shapefiles (or other spatial files)
extent_list <- list.files(system.file("testdata/lakeshapes/", package = "MODIStsp"), full.names = TRUE, "\\.shp$")  

extent_list

# --> Loop on the list of spatial files and run MODIStsp using each of them to 
# automatically define the output extent (A separate output folder is created for
# each input spatial file).

for (single_shape in extent_list) {
  MODIStsp(
    gui = FALSE, 
    opts_file = opts_file, 
    spatmeth = "file", 
    spafile = single_shape, 
    verbose = FALSE, 
    parallel = FALSE
  )
}

# output files are placed in separate folders: 

outfiles_garda <- list.files(
  file.path(tempdir(), "MODIStsp/garda_lake/VI_16Days_1Km_v61/EVI"),
  full.names = TRUE)
outfiles_garda

library(raster)       
plot(raster(outfiles_garda[1]))

outfiles_iseo <- list.files(
  file.path(tempdir(), "MODIStsp/iseo_lake/VI_16Days_1Km_v61/EVI"),
  full.names = TRUE)
outfiles_iseo

plot(raster(outfiles_iseo[1]))

Scheduled Processing

Stand-alone non-interactive execution can be used to periodically and automatically update the time series of a selected product over a given study area. To do that, you should simply:

1. Open the {MODIStsp} GUI, define the parameters of the processing specifying a date in the future as the "Ending Date" and save the processing options. Then quit the program.

2. Schedule non-interactive execution of the launcher installed as seen before (or located in the subdirectory "MODIStsp/ExtData/Launcher" of your library path) as windows scheduled task (or linux "cron" job) according to a specified time schedule, specifying the path of a previously saved Options file as additional argument.

On Linux

1. Edit your crontab by opening a terminal and type:

  crontab -e

1. Add an entry for the launcher. For example, if you have installed it in /usr/bin and you want to run the tool every day at 23.00, add the following row:

  0 23 * * * /bin/bash /usr/bin/MODIStsp -g -s "/yourpath/youroptions.json"

On Windows

1. Create a Task following these instructions; add the path of the MODIStsp.bat launcher as Action (point 6), and specifying -g -s "X:/yourpath/youroptions.json" as argument.

Outputs Format and Naming Conventions

Single-band outputs

Output raster files are saved in specific subfolders of the main output folder. In particular, a separate subfolder is created for each processed original MODIS layer, Quality Indicator or Spectral Index. Each subfolder contains one image for each processed date, created according to the following naming conventions:

myoutfolder/"Layer"/"ProdCode"_"Layer"_"YYYY"_"DOY"."ext"

(e.g., myoutfolder/NDVI/MOD13Q1_NDVI_2000_065.dat)

, where:

• Layer is a short name describing the dataset (e.g., b1_Red, NDII, UI);
• ProdCode is the code name of the MODIS product from which the image was derived (e.g., MOD13Q1);
• YYYY and DOY correspond to the year and DOY (Day of the Year) of acquisition of the original MODIS image;
• ext is the file extension (.tif for GTiff outputs, or .dat for ENVI outputs).

---

Virtual multi-band outputs

ENVI and/or GDAL virtual time series files and RasterStack RData objects are instead stored in the "Time_Series" subfolder if required.

Naming convention for these files is as follow:

<path_of_out_folder>/Time_Series/<vrt_type>/<Sensor>/<Layer>/<ProdCode>_<Layer>_<StartDOY>_<StartYear>_<EndDOY>_<EndYear>_<suffix>.<ext> 

(e.g., /my/out/folder/Time_Series/RData/Terra/NDVI/MOD13Q1_MYD13Q1_NDVI_49_2000_353_2015_RData.RData)

, where:

• <vrt_type> indicates the type of virtual file ("RData", "GDAL" or "ENVI_META");
• <Sensor> indicates to which MODIS sensor the time series belongs ("Terra", "Aqua", "Mixed" or "Combined" (for MCD* products));
• <Layer> is a short name describing the dataset (e.g., "b1_Red", "NDII", "UI");
• <ProdCode> is the code name of the MODIS product from which the image was derived (e.g., "MOD13Q1");
• <StartDOY>, <StartYear>, <EndDOY> and <EndYear> indicate the temporal extent of the time serie created;
• <suffix> indicates the type of virtual file ("ENVI", "GDAL" or "RData");
• <ext> is the file extension (".vrt" for GDAL virtual files, "META" for ENVI meta files or "RData" for R raster stacks).

Accessing the processed time series from R

Preprocessed MODIS data can be retrieved within R either by accessing the single-date raster files, or by loading the saved RasterStack objects.

Any single-date image can be accessed by simply opening it with a raster command:

```{R eval=FALSE, tidy=TRUE, highlight=TRUE} library(raster) modistsp_file <- "/my_outfolder/EVI/MOD13Q1_2005_137_EVI.tif" my_raster <- raster(modistsp_file)

`rasterStack` time series containing all the processed data for a given parameter (saved in the `"Time Series/RData"` subfolder - see [here](https://docs.ropensci.org/MODIStsp/articles/output.html) for details) can be opened by: 

```r
in_virtual_file <- "/my_outfolder/Time_Series/RData/Terra/EVI/MOD13Q1_MYD13Q1_NDVI_49_2000_353_2015_RData.RData" 
indata          <- get(load(in_virtual_file))

This second option allows accessing the complete data stack and analyzing it using the functionalities for raster/raster time series analysis, extraction and plotting provided for example by the {raster} or {rasterVis} packages.

Extracting Time Series Data on Areas of Interest

{MODIStsp} provides an efficient function (MODIStsp_extract()) for extracting time series data at specific locations. The function takes as input a RasterStack virtual object created by MODIStsp() (see above), the starting and ending dates for the extraction and a standard Sp* object (or an ESRI shapefile name) specifying the locations (points, lines or polygons) of interest, and provides as output a xts object or data.frame containing time series data for those locations.

If the input is of class SpatialPoints, the output object contains one column for each point specified, and one row for each date. If it is of class SpatialPolygons (or SpatialLines), it contains one column for each polygon (or each line), with values obtained applying the function specified as the FUN argument (e.g., mean, standard deviation, etc.) on pixels belonging to the polygon (or touched by the line), and one row for each date.

As an example the following code:

  #Set the input paths to raster and shape file
  infile    <- 'myoutfolder/Time_Series/RData/Mixed/MOD13Q1_MYD13Q1_NDVI_49_2000_353_2015_RData.RData'  
  shpname   <- 'path_to_file/rois.shp'  
  #Set the start/end dates for extraction
  startdate <- as.Date("2010-01-01")  
  enddate   <- as.Date("2014-12-31")    
  #Load the RasterStack
  inrts     <- get(load(infile)) 
  # Compute average and St.dev
  dataavg   <- MODIStsp_extract(inrts, shpname, startdate, enddate, FUN = 'mean', na.rm = T)
  datasd    <- MODIStsp_extract (inrts, shpname, startdate, enddate, FUN = 'sd', na.rm = T)
  # Plot average time series for the polygons
  plot.xts(dataavg) 

loads a RasterStack object containing 8-days 250 m resolution time series for the 2000-2015 period and extracts time series of average and standard deviation values over the different polygons of a user's selected shapefile on the 2010-2014 period.

Problems and Issues

Solutions to some common installation and processing problems can be found in {MODIStsp} FAQ: https://docs.ropensci.org/MODIStsp/articles/faq.html

Please report any issues you may encounter in our issues page on GitHub: https://github.com/ropensci/MODIStsp/issues

Citation

To cite MODIStsp please use:

L. Busetto, L. Ranghetti (2016) MODIStsp: An R package for automatic preprocessing of MODIS Land Products time series, Computers & Geosciences, Volume 97, Pages 40-48, ISSN 0098-3004, https://doi.org/10.1016/j.cageo.2016.08.020, URL: https://github.com/ropensci/MODIStsp.

References

lbusett/MODIStsp documentation built on Oct. 16, 2023, 6:59 a.m.