if (FALSE) {
# Extract average and standard deviation values from a rts object created by
# MODIStsp for each polygon of a shapefile, for each date in the period
# between 2001-01-01 and 2014-12-31
# The example uses tif files in testdata/VI_16Days_500m_v61 to build
# a MODIStsp rasterStack corresponding to the 2016 time series of the NDVI index
# over the Como Lake (Italy). It then extracts data on polygons corresponding
# to different land cover classes saved in testdata/extract_polys.shp
# First, prepare the test dataset.
# __NOTE__ To avoid redownloading, here we copy some test data from MODIStsp
# installation folder to tempdir and use it to create a test time series.
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"))
opts_file <- system.file("testdata/test_extract.json", package = "MODIStsp")
MODIStsp(opts_file = opts_file, gui = FALSE, verbose = FALSE)
# Now load the MODIStsp stack: This is a MODIS NDVI time series ranging between
# 2016-01-01 and 2016-12-18
# __NOTE__: MODIStsp rasterStack files are always saved in the "Time_Series\/RData"
# subfolder of your main output folder - see
# "https://docs.ropensci.org/MODIStsp/articles/output.html")
# Specify the filename of the RData RasterStack of interest
stack_file <- file.path(tempdir(),
"MODIStsp/VI_16Days_500m_v61/Time_Series/RData/Terra/NDVI",
"MOD13A1_NDVI_1_2016_353_2016_RData.RData")
basename(stack_file)
ts_data <- get(load(stack_file))
ts_data
# Now load a shapefile containing polygons from which we want to extract data
polygons <- sf::st_read(system.file("testdata/extract_polys.shp",
package = "MODIStsp"), quiet = TRUE)
polygons
# Finally, extract the average values for each polygon and date and plot the
# results
out_dataavg <- suppressMessages(MODIStsp_extract(ts_data, polygons, id_field = "lc_type",
small = FALSE))
head(out_dataavg)
plot(out_dataavg, legend.loc = "topleft")
# use a different summarization function
out_datasd <- MODIStsp_extract(ts_data, polygons, id_field = "lc_type",
FUN = "sd", small = FALSE)
head(out_datasd)
# (See also https://docs.ropensci.org/MODIStsp/articles/Analyze.html for a
# worked-out example)
}
Run the code above in your browser using DataLab