0 
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.

Improve this question
2
  • 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) Commented 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? Commented Apr 3, 2019 at 8:36

1 Answer 1

Reset to default
0

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)
Improve this answer
1
  • 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 Commented Jul 7, 2017 at 19:37

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.