I'm currently running Thornthwaite Evapotranspiration calculation for huge raster stacks. I'm following a solution by Robert Hijmans from another stack overflow thread: Thornthwaite evapotranspiration on a raster dataset . Error formula not vectorised
The short reproducible example from said thread is as follows:
library(SPEI)
library(raster)
library(zoo)
tm = array(20,c(3,4,12*64))
tm = brick(tm)
dates=seq(as.Date("1950-01-01"), as.Date("2013-12-31"), by="month")
tm<- setZ(tm,dates)
names(tm) <- as.yearmon(getZ(tm))
#thornthwaite ET
th <- function(Tave, lat) {
as.vector(SPEI::thornthwaite(as.vector(Tave), lat))
}
a <- raster(tm)
lat <- init(a, "y")
#out <- raster::overlay(tm, lat, fun = th)
out <- brick(tm, values=FALSE)
for (i in 1:ncell(tm)) {
out[i] <- th(tm[i], lat[i])
}
The problem, however, is that the cell-wise calculation of Thornthwaite PET seems rather inefficient, and therefore I strive for a solution using overlay or calc which I assume promises more computational efficiency. In the aforementioned stackoverflow thread another solution using overlay was introduced - however, I fail to reproduce it.
Any solution potentially increasing computational efficiency would be very welcome!
