Spatial Data Analysis: Introduction to Raster Processing: Part 4

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! 😄 👍 👍