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.


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.

1 Answer 1


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] ---


# 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)
  • Thanks, your code works on the test raster stack (except I removed the second comma after scale = 2). The first layer in the result is NA because the two month lag does not work for the first month in the series, so that makes sense. But when I try this with my data I get the same error as above: Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : cannot use this function.
    – jk22
    Mar 28, 2018 at 13:45
  • Am I getting this error because calc() works on chunks of the data but actually needs to have access to the entire series? If so, how can I do this without calc()? The raster stack I'm trying to use has dimensions 585 x 1386 x 360
    – jk22
    Mar 28, 2018 at 14:15
  • Hi! Actually calc works on the entire series in a per pixel basis (from the help you get this: "If x is a RasterStack or RasterBrick, fun should operate on a vector of values (one vector for each cell)"). The problem here is that the function seems not to be working on your case (which is difficult to check what is happening exactly without accessing the data). Try this out: funSPEI <- function(x, scale=2, na.rm=TRUE,...) as.numeric((spei(x, scale=scale, na.rm=na.rm, ...))$fitted); rstSPEI <- calc(rstack, fun = funSPEI)
    – Kamo
    Mar 28, 2018 at 18:02
  • Hi again! I was able to generate your error and it is happening because of NA's in the series so, if you turn the na.rm=TRUE, the problem seems to be fixed. Check out the edited solution above.
    – Kamo
    Mar 28, 2018 at 18:26
  • I have run into calc() and overlay() failing during the internal check, for no apparent reason, so many times that, if possible, I coerce into my stack into a SpatialPixelsDataFrame object and then use apply() to operate on rows. Mar 28, 2018 at 19:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.