Торнвэйтское суммарное испарение в наборе растровых данных.Формула ошибки не векторизована

Question

Я пытаюсь рассчитать суммарное испарение (ET) для запуска SPEI для набора растровых данных, используя формулу Thornthwaite ET, входящую в пакет SPEI

, это мой код

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

tm = array(1:(3*4*12*64),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) {
  SPEI::thornthwaite(Tave, lat)
} 

lat <- setValues(a, coordinates(tm)[, "y"])
out <- raster::overlay(tm, lat, fun = th)

, ноя получил следующую ошибку:

Error in (function (x, fun, filename = "", recycle = TRUE, forcefun = FALSE,  : 
  cannot use this formula, probably because it is not vectorized

Пожалуйста, вы можете помочь?

Спасибо миллион

Answer 1

Чтобы сделать эту работу в растре :: оверлей, требуется еще один шаг.Это «векторизация» функции, как было предложено в сообщении об ошибке.Vectorize - это базовая функция, которая создает оболочку для данной функции для работы с mapply.

Пожалуйста, смотрите код ниже:

library(raster)
library(zoo)

tm = array(1:(3*4*12*64),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))
} 

lat <- init(raster(tm), "y")

## now run overlay with Vectorize in the function call
out <- overlay(tm, lat, fun = Vectorize(th))

plot(out)

Вот график из моего тестового прогона

Answer 2

Я не совсем уверен, почему это не удается.Вот обходной путь

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])
}