### 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`

or`RasterBrick`

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)
```

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