[GIS] Using spei() function on time series from rasterstack in R

rrastertime series

I want to use the spei() function in the SPEI package on a raster stack of monthly time series of water balance data. The stack has 30 years of monthly data, so 360 layers. The documentation for spei() https://cran.r-project.org/web/packages/SPEI/SPEI.pdf states the input should be "a vector, matrix or data frame with time ordered values", along with a lag (number of months). Can I run spei() without converting each pixel in the stack to a time series vector? I have explored this thread -> Apply SPEI::hargreaves function to time series from each pixel rasterbrick R which got me started setting dates as the Z dimension in the stack. 

library(raster)
library(SPEI)
library(zoo)

wb.files <- list.files(path = <pathway to files>, full.names = TRUE)
wb.stack <- stack(wb.files)
dates <- seq(as.Date("1971-01-01"), as.Date("2000-12-31"), by="month")
wb.stack <- setZ(wb.stack, dates)
names(wb.stack) <- as.yearmon(getZ(wb.stack))

Then when I try to run the function, I get an error.

system.time(spei.r <- spei(wb.stack, 2))

The error I get is:

#Error in sum(is.na(data)) > 0 & na.rm == FALSE : 
#  operations are possible only for numeric, logical or complex types

In the thread above, the raster::overlay() function is used, and the equivalent for a single raster stack is calc() so I tried using that too.

speicalc <- function(dat, sc) {
  SPEI::spei(dat, sc)
} 

system.time(spei.r <- raster::calc(wb.stack, 2, fun = speicalc))

The error I get in that case is:

#Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : cannot use this function

I would like to avoid collapsing every pixel to a vector or data frame, which seems like it would be very inefficient.

Best Answer

The output of spei() function is an object of class spei with multiple elements and as such, calc() does not know what to do with it.

So, guessing that you want to extract the $fitted part of the output of spei [i.e., fitted: time series with the values of the Standardized Precipitation-Evapotranspiration Index (SPEI) or the Standardized Precipitation Index (SPI)] you need to extract only that part (which is a ts object) convert it to a numeric vector that calc likes and also pass the scale parameter in the modified spei function. It gets like this:

--- [Edited in 28/03/2018] ---

library(SPEI)
library(raster)

# Generate a sample raster stack time series with 360 layers (this will take a while...)
r <- raster(nrows=10,ncols=10,vals=rnorm(100))
rstack <- stack(r)
for(i in 1:359){
  rstack <- stack(rstack, raster(nrows=10,ncols=10,vals=rnorm(100)))
  cat(paste("..",round(((i+1)/360)*100,1),"%")) # check progress in %
}

# Change the function SPEI so it outputs a numeric vector from fitted
# Pass scale value
rstSPEI <- calc(rstack, fun = function(x, scale = 2, na.rm = TRUE, ...) as.numeric((spei(x, scale = scale, na.rm = na.rm, ...))$fitted))

# Or using your example
spei.r <- raster::calc(wb.stack, fun = function(x, scale = 2, na.rm = TRUE, ...) as.numeric((spei(x, scale = scale, na.rm = na.rm, ...))$fitted))

[Note: I am not sure why but, the first image returned by spei in the example is NA... (explained!: check comments below)]

You can also declare the function outside and use it in calc making a specific reference to the na.rm parameter (which usually gives some problems) like this:

# Declare the function to use
funSPEI <- function(x, scale=2, na.rm=TRUE,...) as.numeric((spei(x, scale=scale, na.rm=na.rm, ...))$fitted)

rstSPEI <- calc(rstack, fun = funSPEI)
Related Question