# Applying SPEI::hargreaves function to time series from each pixel rasterbrick using R?

library(SPEI)
library(raster)


Basically, given a monthly time series for two variables a and b as follows:

a = array(1:(3*4*12*64),c(3,4,12*64))
a = brick(a)

dates=seq(as.Date("1950-01-01"), as.Date("2013-12-31"), by="month")
a<- setZ(a,dates)
names(a) <- as.yearmon(getZ(a))

b = array(1:(3*4*12*64),c(3,4,12*64))
b = brick(b)

dates=seq(as.Date("1950-01-01"), as.Date("2013-12-31"), by="month")
b<- setZ(b,dates)
names(b) <- as.yearmon(getZ(b))


Both a and b are time series with a time dimension. Now, I would like to apply the function SPEI::hargreaves to time series of each pixel in a and b and return a rasterbrick C

From the SPEI package , har <- hargreaves(TMIN,TMAX,lat=37.6475). Here, consider Tmin=a, Tmax=b and lat=latitude of each pixel which is the same in a and b. I will parallelize the process once I get an idea how to apply the function to my rasterbricks.

At the moment, I am using data.table to collapse my rasterbricks and then applying hargreaves to the table. This approach is very inefficient so far.

• The solution above works for SPEI::hargreaves but not for SPEI::thornthwaite, which is an alternative formula to calculate PET. The solution involves using: thor <- Vectorize(SPEI::thornthwaite) Instead of the function(){...} wrapper (@sboysel, I'm not sure what that was for anyway? would be really helpful if you could explain) Aug 14, 2018 at 20:13
• I have exactly the same problem and tried the code of sboysel but I get an error because the formula is not vectorized. Anyone an idea how to solve this?
– GVDB
Apr 3, 2019 at 8:36

See ?raster::overlay. You could try defining a new brick containing the latitudes of each cell:

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

a = array(1:(3*4*12*64),c(3,4,12*64))
a = brick(a)

dates=seq(as.Date("1950-01-01"), as.Date("2013-12-31"), by="month")
a<- setZ(a,dates)
names(a) <- as.yearmon(getZ(a))

b = array(1:(3*4*12*64),c(3,4,12*64))
b = brick(b)

dates=seq(as.Date("1950-01-01"), as.Date("2013-12-31"), by="month")
b<- setZ(b,dates)
names(b) <- as.yearmon(getZ(b))

har <- function(Tmin, Tmax, lat) {
SPEI::hargreaves(Tmin, Tmax, lat)
}

lat <- setValues(a, coordinates(a)[, "y"])
c <- raster::overlay(a, b, lat, fun = har)

• what could be the reason behind tgis warning message after I run your code above: Warning message: In if (class(test2) == "try-error" | length(test2) < 5) { : the condition has length > 1 and only the first element will be used Jul 7, 2017 at 19:37