Background
In the fourth part of this tutorial series on Spatial Data Analysis using the raster
package, we will explore more functionalities, this time related to time-series analysis of raster data. For more information on raster data processing, see here, as well as the tutorial part-1, tutorial part-2, and, tutorial part-3, of this series.
We will use an Enhanced Vegetation Index (EVI), 5-year time-series (from the year 2012 to 2016) from Terra/MODIS satellite/sensor platform for the Peneda-Geres National Park (PGNP, in NW Portugal) to develop some examples.
This data corresponds to MODIS’s MOD13Q1 data product version-006, (+info here) which has 250m of spatial resolution and 16-days of temporal resolution (more precisely, the product is generated by maximum daily value composites for each 16-day period). This means that each year has a total of 23 observations. This data was downloaded from the EarthData platform and later assembled and re-projected to WGS 1984 – UTM 29N Coordinate Reference System (CRS) using the MODIS Reprojection Tool – MRT (sorry, but these pre-processing steps are outside the scope of this tutorial 😁).
In this post, we will introduce RasterBrick
‘s, a multi-layer raster object typically created from a multi-layer (or multi-band) file, although they can also exist entirely in memory. These objects are similar to RasterStack
s, but processing time should be shorter when using a RasterBrick
(irrespective, if values are on disk or in memory). However, these objects are less flexible, as they can only point to a single file, while RasterStacks
can point to multiple files.
Besides the raster
package, we will also work with rts
, which provides classes and methods for manipulating and processing raster time-series data (e.g. a time-series of satellite images). A raster time-series object is created by combining a RasterStack
or RasterBrick
object (from raster
package) and a set of dates of class POSIXct
, POSIXt
, Date
, timeDate
. The time information in rts
is then handled by a xts
object.
The function rts
is used to build a raster time-series (either a RasterBrickTS
or RasterStackTS
) which is simply an S4 object composed of two slots:
- Slot [raster]: a
RasterStack
orRasterBrick
object. - Slot [time]: a
xts
object with dates for each layer in the raster object.
One key advantage of using the rts
package is to facilitate the subset, extraction, or application of functions by specific periods using date notation (instead of integer or name indices like in raster
).
Making a Raster Time-series Object With RTS Package
First up: download and uncompress the sample data! 👍 The .zip archive contains a multi-layer GeoTIFF file with a 16-day (composite) EVI time-series from 2012 to 2016. This means that we have a total of 23 images per year and a total of 115 layers in the file (for the whole five-year series).
## Create a folder named data-raw inside the working directory to place downloaded data
if(!dir.exists("./data-raw")) dir.create("./data-raw")
## If you run into download problems try changing: method = "wget"
download.file("https://raw.githubusercontent.com/joaofgoncalves/R_exercises_raster_tutorial/master/data/MODIS_EVI_TS_PGNP_MultiBand.zip", "./data-raw/MODIS_EVI_TS_PGNP_MultiBand.zip", method = "auto")
## Uncompress the zip file
unzip("./data-raw/MODIS_EVI_TS_PGNP_MultiBand.zip", exdir = "./data-raw")
Creating the RasterBrick
from the downloaded data is really easy and similar to using the stack
. The main advantage is that the input is just one multi-layer file instead of a vector with multiple files as usual when using stack
. We will also change the names of each layer for making them more readable. Let’s see how this goes:
library(raster)
# Load the raster data into a RasterBrick object
rst <- brick("./data-raw/MOD13Q1.2012_2016.PGNP_250m_EVI_16days.tif")
names(rst) <- paste("EVI",1:nlayers(rst),sep="_")
MODIS data products, such as MOD13Q1 (and others alike), use 7-digit long dates composed by the year, (4 digits) followed by the Julian day (3 digits) to identify the reference date of an image composite. For example, the date code 2012001 corresponds to 2012-01-01 (in YYYY-mm-dd format), and, 2012161 to 2012-06-09. Usually, these dates are inscribed in image files or meta-data but, since we don’t have them, we will generate them first and then process them so we can obtain a Date
object for each layer. We can then use these properly formatted dates in the rts
function to create a raster time-series (see ?rts
for more details).
padIt <- function(x)
if(x<10) paste("00",x,sep="") else if(x<100 & x>=10) paste("0",x, sep="") else return(as.character(x))
padWithZeros <- function(x)
sapply(x, FUN = padIt)
# Generate the MODIS-like dates for each layer
julDay <- padWithZeros(rep(seq(from = 1,to = 365,by = 16),5))
yrs <- as.character(rep(2012:2016, each = 23))
MODISYrJday <- paste(yrs, julDay, sep="")
# Print out the MODIS dates for the year 2012
print(MODISYrJday[1:23])
## [1] "2012001" "2012017" "2012033" "2012049" "2012065" "2012081" "2012097"
## [8] "2012113" "2012129" "2012145" "2012161" "2012177" "2012193" "2012209"
## [15] "2012225" "2012241" "2012257" "2012273" "2012289" "2012305" "2012321"
## [22] "2012337" "2012353"
Now that we have our MODIS-like dates (in year and Julian day format) we need to convert them into a more ‘human-readable’ format, also accepted by rts
function:
# Extract the year
MOD.getYear<-function(x)
as.integer(sapply(x,FUN=function(x) substr(x,1,4)))
# Extract de julian day
MOD.getDOY<-function(x)
as.integer(sapply(x,FUN=function(x) substr(x,5,7)))
# Process the MODIS-like date into YYYY-mm-dd format as Date object
MOD.getDate<-function(x)
as.Date(sapply(x,FUN=function(x) as.character(as.Date(MOD.getDOY(x)-1,origin=paste(MOD.getYear(x),"01-01",sep="-")))))
MODdates <- MOD.getDate(MODISYrJday)
class(MODdates)
## [1] "Date"
# Check the result for year 2012
print(MODdates[1:23])
## [1] "2012-01-01" "2012-01-17" "2012-02-02" "2012-02-18" "2012-03-05"
## [6] "2012-03-21" "2012-04-06" "2012-04-22" "2012-05-08" "2012-05-24"
## [11] "2012-06-09" "2012-06-25" "2012-07-11" "2012-07-27" "2012-08-12"
## [16] "2012-08-28" "2012-09-13" "2012-09-29" "2012-10-15" "2012-10-31"
## [21] "2012-11-16" "2012-12-02" "2012-12-18"
With the date vector (class Date
) for each layer, we can now generate a raster time-series with rts
constructor:
# Install the rts package
if(!("rts" %in% installed.packages()[,1]))
install.packages(c("rts"), dependencies = TRUE)
library(rts)
rstTS <- rts(rst, MODdates)
Sub-setting Raster Time-series
With the RasterBrickTS
object created, we can extract sub-sets of the data for particular periods or dates (check ?rts::subset
for more details):
# Subset a specific period
rstTSsubset1 <- subset(rstTS,"2013-05-15/2014-08-25")
# Subset the whole year of 2012
rstTSsubset2 <- subset(rstTS,"2012")
# Subset years from 2013 to 2014
rstTSsubset3 <- subset(rstTS,"2013/2014")
# Subset all years from (and including) 2014 to the series end
rstTSsubset4 <- subset(rstTS,"2014/")
# Subset all to the end of 2014
rstTSsubset5 <- subset(rstTS,"/2014")
# Subset all to the end of July 2014
rstTSsubset6 <- subset(rstTS,"/2014-06")
# Subset a specific month
rstTSsubset7 <- subset(rstTS,"2016-05")
# Plot the May 2016 data
plot(rstTSsubset7)
[Intermediate] Spatial Data Analysis with R, QGIS & More. this course you will learn how to:
- Work with Spatial data and maps
- Learn about different tools to develop spatial data next to R
- And much more

As you can see, the subset
function is pretty handy for extracting parts of a time-series. The date parameter format must be left-specified with respect to the standard ISO:8601 time format “CCYY-MM-DD HH:MM:SS.” It is also possible to specify a range of times via the index-based sub-setting, using ISO-recommended “/” as the range operator. Generally, it works with “from/to” dates, where using both is optional. If one side is missing, it is interpreted as a request to retrieve layers from the beginning or through the end of the raster time series.
Apply Functions Over the Time-series
One of the best applications of the rts
package for processing raster time-series is the ability to apply a specified function to distinct periods. Let’s see how we can do this:
# Apply function to each quarter
# Mean
rstTS_quarterlyMN <- apply.quarterly(rstTS, FUN = mean, na.rm=TRUE)
# Standard-deviation
rstTS_quarterlySD <- apply.quarterly(rstTS, FUN = sd, na.rm=TRUE)
# Apply function to each year
# Mean
rstTS_yearlyMN <- apply.yearly(rstTS, FUN = mean, na.rm=TRUE)
# Standard-deviation
rstTS_yearlySD <- apply.yearly(rstTS, FUN = sd, na.rm=TRUE)
# Plot the time-series for annual EVI
plot(rstTS_yearlyMN)

As we can see, these functions can be very useful for applying specific functions over certain calendar periods without the hassle of having to specify indices – you simply work with dates, which are nicer! 👍 😉
However, in certain cases, you may want to work with the “simpler” / more general functions that the raster
package offers to apply functions over a raster time-series. In that case, calc
and stackApply
can be used.
The difference between these functions is that calc
applies the defined function over the whole series pixel-by-pixel, while stack
applies a function on sub-sets of a RasterStack
or RasterBrick
. For this function, the layers to be combined are indicated by an integer vector with indices. The function used should return a single value, and the number of layers in the output Raster*
equals the number of unique values in indices. In the opposite hand, calc
allows having a function that outputs multiple values and, in that case, a multi-layer Raster*
object is returned with one layer per output value.
Also, keep in mind that for large objects, calc
will compute values chunk by chunk. This means that for the result of fun to be correct, it should not depend on having access to all values at once. Let’s grab some examples.
Using the calc
function to provide a global average and standard-deviation of the entire raster time-series:
rstMean <- calc(rst, fun = mean)
rstStd <- calc(rst, fun = sd)
Now, let’s use calc
for a multi-value output with a specific function:
# Calculate quantiles
# NOTE: you have to use na.rm=TRUE to make this work
rstQuantiles <- calc(rst, fun = function(x,...) as.numeric(quantile(x,probs=c(0.05, 0.5, 0.95),...)), na.rm=TRUE)
print(rstQuantiles)
## class : RasterBrick
## dimensions : 186, 179, 33294, 3 (nrow, ncol, ncell, nlayers)
## resolution : 250, 250 (x, y)
## extent : 549486, 594236, 4613206, 4659706 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer.1, layer.2, layer.3
## min values : -290.7, 13.0, 271.5
## max values : 4229.3, 5572.0, 7458.0
Et voilà! ❕ We have three layers; one for each calculated quantile.
Now, switching to stackApply
, we will emulate the behavior of rts apply.yearly
function:
rstYrMean <- stackApply(rst, fun=mean, indices = rep(1:5,each=23))
Well… To be honest, I was not aware of this but, strangely the rts
function took much more time to calculate the annual averages… Check out the comparison below:
# stackApply with RasterBrick
system.time({rstYrMean <- stackApply(rst, fun=mean, indices = rep(1:5,each=23))})
## user system elapsed
## 0.52 0.04 0.58
# apply.yearly with RasterBrickTS
system.time({rstTS_yearlyMN <- apply.yearly(rstTS, FUN = mean, na.rm=TRUE)})
## user system elapsed
## 23.47 1.25 24.75
Let’s check if the data is equal to be sure:
for(i in 1:nlayers(rstYrMean))
print(compareRaster(rstYrMean[[i]], rstTS_yearlyMN@raster[[i]], values=TRUE))
## [1] TRUE
## [1] TRUE
## [1] TRUE
## [1] TRUE
## [1] TRUE
Yup, all the layers are equal…
Not sure why this is happening though… 🤔. Do you have some ideas/comments on this?
This concludes our exploration of the raster and the rts packages for this post. We have covered a couple of useful things for processing raster time-series. Hope you find it useful! 😄 👍 👍
Leave a Reply