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

knitr::opts_chunk$set(echo = TRUE)
library(rmarkdown)
library(knitr)
library(raster)
library(base)
library(sf)
library(MODIStsp)
library(dplyr)
test_zip <-  system.file("testdata/VI_16Days_500m_v6/NDVI.zip",
                         package = "MODIStsp")
dir.create(file.path(tempdir(), "MODIStsp/VI_16Days_500m_v61"),
           showWarnings = FALSE, recursive = TRUE)
utils::unzip(test_zip, 
             exdir = file.path(tempdir(), "MODIStsp/VI_16Days_500m_v61"))
Sys.setlocale("LC_TIME", "en_US.UTF-8")

Accessing the raster time series

Preprocessed MODIS data can be retrieved within R either by accessing the single-date raster files, or by loading the saved RasterStack objects (see here for a description of {MODIStsp} output folder structure and naming conventions).

To test this functionality, you can run the following example:

library(MODIStsp)
opts_file <- system.file("testdata/test_extract.json", package = "MODIStsp")
MODIStsp(opts_file = opts_file, gui = FALSE, verbose = FALSE)

This will download a yearly time series of MODIS NDVI data and subset it over the region of the Como Lake in Italy.

After the download and processing finishes (it will take a while, depending on your network speed), the MODIS time series will be placed in subfolder MODIStsp/VI_16Days_500m_v61 of R tempdir(). In this case, I have them in "file.path(tempdir(), "VI_16Days_500m_v61")". If you want to save them elsewhere, run MODIStsp(opts_file = opts_file, gui = TRUE) and select a different output folder.

Single-date processed rasters can be accessed by simply opening them with a raster command:

```{R echo=TRUE, message=FALSE, warning=FALSE, highlight=TRUE, tidy=TRUE} modistsp_file <- file.path(tempdir(),"MODIStsp/VI_16Days_500m_v61/NDVI", "MOD13A1_NDVI_2016_177.tif") modistsp_file

my_raster <- raster(modistsp_file) plot(my_raster)

`RasterStack` (or GDAL vrt) time series containing all the processed data for a
given parameter are saved in the `Time Series/RData` subfolder, and can be accessed
by: 

```r
# Load the NDVI time series
ts_folder <- file.path(tempdir(), "MODIStsp/VI_16Days_500m_v61/Time_Series")
in_virtual_file <- file.path(ts_folder, "RData/Terra/NDVI/MOD13A1_NDVI_1_2016_353_2016_RData.RData")
in_virtual_file

ts_data   <- get(load(in_virtual_file))
ts_data

, which gives us a 23-band RasterStack in the ts_data variable.

This RasterStack can be analyzed using the functionalities for raster/raster time series analysis, extraction and plotting provided for example by the {raster} or {rasterVis} packages:

# plot some dates
plot(ts_data[[c(1,5,10,15,20)]], axes = FALSE, horizontal = T)

# Extract one date from the stack

mydate   <- as.Date("2016-01-01")
substack <- subset(ts_data, which(getZ(ts_data) == mydate)) %>% setZ(mydate)
substack  

# Extract multiple dates from the stack

mindate   <- as.Date("2016-01-01")
maxdate   <- as.Date("2016-04-01")
substack  <- subset(ts_data, 
                    which(getZ(ts_data) >= mindate & getZ(ts_data) <= maxdate))
substack  

# Compute monthly averages

month_avg <- stackApply(ts_data, months(getZ(ts_data)), fun = mean)
month_avg

plot(month_avg, axes = FALSE, horizontal = T)

---

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.

To test MODIStsp_extract() we can start by loading the NDVI RasterStack time series created in the previous example:

ts_folder <- file.path(tempdir(), "MODIStsp/VI_16Days_500m_v61/Time_Series")
in_virtual_file <- file.path(ts_folder, "RData/Terra/NDVI/MOD13A1_NDVI_1_2016_353_2016_RData.RData")
ts_data   <- get(load(in_virtual_file))
ts_data

, which again gives us the 23-band RasterStack in the ts_data variable.

To extract the NDVI data over selected areas, we can then use the MODIStsp_extract() function as follows:

# load a shapefile containing polygons --> Here we use a test shapefile in 
# the testdata folder of MODIStsp.

polygons <- sf::st_read(system.file("testdata/extract_polys.shp",
                                       package = "MODIStsp"),
                           quiet = TRUE)
polygons

# Now extract the average values for each polygon and date and plot the results

out_dataavg <- MODIStsp_extract(ts_data, polygons, id_field = "lc_type", 
                                small = FALSE)

head(out_dataavg)

# Other summarization functions can be used, by specifying the "FUN" argument. 
# Compute the Standard Deviation over each polygon: 
out_datasd <- MODIStsp_extract(ts_data, polygons, id_field = "lc_type",
                               FUN = "sd", small = FALSE)

The output is a xts object, with one column for each polygon of the input shapefile, and one row for each date. We can easily plot the computed averages using:

library(xts)
plot(out_dataavg, legend.loc = "topleft")

We can also transform the output dataset to a long format and use {ggplot2} for plotting and other tidyverse-friendly functions for analysis:

library(ggplot2)
out_dataavg_df <- data.frame(date = index(out_dataavg), coredata(out_dataavg)) %>% 
  tidyr::gather(Polygon, Value, -date)

ggplot2::ggplot(out_dataavg_df, aes(x = date, y = Value, color = Polygon)) + 
  geom_line() + theme_light()

require(tibbletime)
require(dplyr)

# Compute monthly averages using tibbletime: 

out_dataavg_tt <- tibbletime::as_tbl_time(out_dataavg_df, index = date) 

month_avg <- out_dataavg_tt %>% 
  dplyr::group_by(Polygon) %>% 
  tibbletime::as_period(period = "monthly")

ggplot2::ggplot(month_avg, aes(x = date, y = Value, color = Polygon)) + 
  geom_line() + theme_light()

Important note: MODIStsp_extract() is usually much faster than the standard raster::extract() function, but does not deal with overlapping polygons. If your polygons are overlapping, please use raster::extract() instead.

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