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

Question

我正在使用作者随可用研究发布的 Rcode 此处计算本研究/科学出版物

的最大气候缺水量

这是一个非常简单的代码，其中输入是我从 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])  
}

矢量化版本（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

Answer 2

我使用 calc() 对自定义函数的实现进行矢量化时遇到了类似的问题，该函数使用 RasterStack 输入上趋势包中的 sens.slope(x) 。这是我的解决方案，希望对其他人有帮助。

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
  cannot use this function

我的输入 RasterStack 包含带有 NA 的单元格，因此我在自定义函数中内置了 NA 处理：

  if(sum(is.na(x)) = length(x)) return(NA)
  if(na.rm) x <- x[!is.na(x)]

但是，如果 x 中的所有值均为 NA，则 sens.slope(x) 会失败，但如果 x 中只有一个非 NA 值，则 sens.slope(x) 也会失败。上面的 NA 处理处理了前一种情况，但不处理后者。更新我的自定义函数以包括 x 中只有一个非 NA 的情况解决了我的问题。

  if(sum(is.na(x)) >= length(x)-1) return(NA)
  if(na.rm) x <- x[!is.na(x)]

这就是我的解决方法，我将几个虚拟数据集推入我的函数中，而不使用 calc() ，直到重现错误：

x <- c(1,2,3,4,5,6) 
myfun(x,na.rm=TRUE)
output was fine

x <- c(1,2,3,4,5,NA) 
myfun(x,na.rm=TRUE) 
output was fine

x <- c(NA,NA,NA,NA,NA,NA) 
myfun(x,na.rm=TRUE) 
output was fine

x <- c(NA,NA,NA,NA,5,NA) 
myfun(x,na.rm=TRUE) 
output was not fine, error