r raster 如何矢量化函数以与叠加层一起使用

Question

我想做一个使用此函数的栅格叠加，其中 x 和 y 是栅格，Tbase 是一个常数。

gddFunction <- function(x, y, z) {(x + y)/2 - z}

喜欢这个

library(raster)
gdd <- overlay(x = tmax_clamped, y = tmin_clamped, z = Tbase, fun = gddFunction)

我收到以下错误消息

(函数 (x, fun, filename = "", recycle = TRUE, forcefun = FALSE, : 不能使用这个公式，可能是因为它不是 矢量化

我尝试按照帮助中的示例进行操作，如下

 vgddFunction <- Vectorize(gddFunction, SIMPLIFY = FALSE)
 gdd <- overlay(x = tmax_clamped, y = tmin_clamped, z = Tbase, fun = vgddFunction)

但收到相同的错误消息。

我怀疑问题是我需要 z 是常数值的向量，或者可能是常数值的栅格层。但我希望 Vectorize 操作能做到这一点。

直截了当的地图代数工作，如

gdd <- (x + y)/2 - z

但它很慢，我希望通过使用叠加来加快计算速度

Answer 1

raster::overlay 仅允许将 raster 参数作为额外参数通过 ... 传递（请参阅 ?overlay）。所以你不能直接做现在正在做的事情。但是我们可以避免这个问题。一个可重现的例子：

library(raster)
r <- raster(ncol = 10, nrow = 10)
r1 <- init(r, fun = runif)
r2 <- init(r, fun = runif)
z <- 3

你的错误：

gddFunction <- function(x, y, z) {(x + y)/2 - z}
gdd <- overlay(x = r1, y = r2, z = z, fun = gddFunction)

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

如果我们不能直接传递z，我们可以通过一些函数式编程间接传递：

gddFunction2 <- function(z) {
  function(x, y) (x + y) / 2 - z
}
gdd <- overlay(x = r1, y = r2, fun = gddFunction2(z))

这里，gddFunction2 将返回另一个函数，该函数只有一个 x 和 y 参数，可以接受 rasters，但您仍然可以更改常量 z，它现在是 @ 的参数987654335@.

你也可以这样写：

gddFunction <- gddFunction2(z)
gdd <- overlay(x = r1, y = r2, fun = gddFunction)

请注意，您原来的 gddFunction 已经矢量化，因此错误消息是在错误地猜测问题。