Use spei() function with user specified parameters from a rasterstack in R?

Question

This is a follow up to Using spei() function on time series from rasterstack in R?

I am trying to use the spei() function from the SPEI package to calculate drought metrics. I have a rasterstack of 30 years of monthly data, so 360 layers. I want to run spei() once on historic climate data to calculate the coefficients for fitting the SPEI index. Then I want to apply those coefficients to new monthly data for recent and future climate to actually calculate SPEI (fitted) values on those data. So, @Kamo helped me in the previous Q&A to run spei() and return a rasterstack.

library(SPEI)
library(raster)

# Generate a sample raster stack time series with 360 layers 
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 of coefficients
# Modified from @Kamo
funSPEI <- function(x, scale=2, na.rm=TRUE) {
  as.numeric((spei(x, scale=scale, na.rm=na.rm))$coefficients)
}

rst.coeff <- calc(rstack, fun = funSPEI)

Now I want to use these coefficients as the param argument in spei(). According to spei documentation (https://cran.r-project.org/web/packages/SPEI/SPEI.pdf), user specified params "should be a three-dimensional array with dimensions (3,number of series in data,12), containing twelve parameter triads (xi, alpha, kappa) for each data series, one for each month".

Instead of 12 parameter triads, what I have is a raster stack of 36 layers (3 parameters times 12 months). How can I use that as input to spei()?

I have tried:

funSPEI.2 <- function (x, scale=2, na.rm=TRUE, p) { 
  as.numeric((spei(x, scale=scale, na.rm=na.rm, p))$fitted) 
}

rstFit <- calc(rstack, fun = funSPEI.2, params  = rst.coeff))

But I get the error:

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

Answer 1

I found another way to do this using a single run of spei() for the entire time series and modifying my function to include the ref.start and ref.end arguments. So multiple calls of spei() are not needed. The technique relies on converting the raster stack to a matrix in which each row is a time series of values. Then, apply is used to run spei() on each row. The resulting matrix is converted then back to a raster stack.

In my case, I'm using the first 30 years of monthly data as the reference period.

@Kamo and @Jeffrey-Evans helped getting me started on this.

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

# Generate a sample raster stack time series with 720 layers
# This would be 60 years of monthly data
r <- raster(nrows=10,ncols=10,vals=rnorm(100))
rstack <- stack(r)
for(i in 1:719){
  rstack <- stack(rstack, raster(nrows=10,ncols=10,vals=rnorm(100)))
  cat(paste("..",round(((i+1)/720)*100,1),"%")) # check progress in %
}

r.mat <- as.matrix(rstack)

# Run spei()
funSPEImat <- function(x, sc, start, end, na.rm=TRUE,...) {
  dat <- ts(x, start = c(1971, 1), end = c(2030, 12), frequency = 12)
  as.numeric((spei(dat, sc, ref.start = start, ref.end = end, na.rm=na.rm, ...))$fitted) 
}

fitted.mat <- t(apply(r.mat, 1, funSPEImat, sc = 2, start = c(1971, 1), 
                                 end = c(2000, 12)))


# Convert back to raster brick
spei <- setValues(rstack, fitted.mat)
dates <- seq(as.Date("1971-01-01"), as.Date("2030-12-31"), by="month")
names(spei) <- as.yearmon(dates)