带有用户定义函数的 calctest(x[1:5]) 中的 raster::calc() 错误

Question

我正在计算此study/scientific publication 的最大气候缺水量，使用作者发布的 Rcode 以及可用的研究here

这是一个非常简单的代码，其中的输入是我从 Terraclimate 产品下载的月降雨量栅格。从每个降雨栅格（月）中减去 100 毫米/月（蒸散量）后，会创建一个堆栈，并在 calc() 中使用以下函数

wd = stack(month.rainfall)-100 # 100 is the evapotranspiration in mm/month

# MCWD Function
mcwd.f = function(x){
result= as.numeric(x)
for(i in 1:length(result)){
wdn = result[i]
wdn1 = result[i-1]

if(i==1){
  if(wdn>0){ result[i]=0}
  else{result[i]=wdn}
}

if(i!=1){
  cwd = wdn1+wdn
  if( cwd < 0){ result[i]=cwd}
  else{result[i]=0}
   }
 }

return(result)  
 }


# Applying the Function
cwd = calc(wd, fun = mcwd.f)

但是，在我运行从 cwd 开始的行之后，我得到了错误 .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) 中的错误： 无法使用此功能

为什么会出现这个错误？我该如何解决？

我查看了类似的帖子-1、2，并使用了 na.rm=TRUE 等，但我仍然遇到相同的错误。我还看到一些帖子说 calc() 不是很好，但我不知道如何解决这个问题或找到替代方案，因为我正在使用我正在复制的研究作者发布的代码。

编辑-我使用了一个包含 100 个随机数（没有 NA）的向量来测试 mcwd.f 函数并且它可以工作。所以我认为这与拥有 NA 像素有关，如果我在 calc() 中使用 na.rm=TRUE，我认为应该可以解决这个问题，但它并没有解决它。

Answer 1

当您提出 R 问题时，请始终包含一个最小的、独立的、可重现的示例。您说您的功能“有效”，但您没有提供输入和预期的输出数据。

我用你的功能得到了这个

mcwd.f(c(-5:5))
# [1]  -5  -9 -12 -14 -15 -15 -14 -12  -9  -5   0

这是你的函数的更简洁的版本

mcwd.f1 <- function(x) {
    x[1] <- min(0, x[1])
    for(i in 2:length(x)){
        x[i] <- min(0, sum(x[c(i-1, i)]))
    }
    return(x)  
}
mcwd.f1(c(-5:5))
# [1]  -5  -9 -12 -14 -15 -15 -14 -12  -9  -5   0

还有一个有助于解释下一个的变体

mcwd.f1a <- function(x) {
    x <- c(0, x)
    for(i in 2:length(x)){
        x[i] <- min(0, sum(x[c(i-1, i)]))
    }
    return(x[-1])  
}

矢量化版本 (after Gabor Grothendieck)

mcwd.f2 = function(x) {
    f <- function(x, y) min(x + y, 0)
    Reduce(f, x, 0, accumulate = TRUE)[-1]
}
mcwd.f2(c(-5:5))
#[1]  -5  -9 -12 -14 -15 -15 -14 -12  -9  -5   0

现在将该函数与 RasterStack 和 calc 一起使用

示例数据：

library(raster)
slogo <- stack(system.file("external/rlogo.grd", package="raster")) 
s <- slogo-100

使用计算：

x <- calc(s, mcwd.f1)
y <- calc(s, mcwd.f2)

这适用于这两个功能（也适用于您的mcwd.f）。尝试使用NA 值

s[1000:3000] <- NA 
x <- calc(s, mcwd.f1)
y <- calc(s, mcwd.f2)

也有效。

x
#class      : RasterBrick 
#dimensions : 77, 101, 7777, 3  (nrow, ncol, ncell, nlayers)
#resolution : 1, 1  (x, y)
#extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
#crs        : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source     : memory
#names      :  red, green, blue 
#min values : -100,  -200, -300 
#max values :    0,     0,    0 

但不适合你的功能

y <- calc(s, mcwd.f)
Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
  cannot use this function

当你这样做时，你就会明白为什么

mcwd.f(c(-5:5,NA))
#Error in if (cwd < 0) { : missing value where TRUE/FALSE needed

为了让你正常工作（但你不应该因为它效率低下），你可以像这样添加第一行

mcwd.f = function(x){
    if (any(is.na(x))) return(x*NA)
    ...

假设当存在一个 NA 时，单元格中的所有值都是或应该成为 NA