# 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) – Andrea Aug 14 '18 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 '19 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` – code123 Jul 7 '17 at 19:37