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")
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)
{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.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.