raster::calc() error in calctest(x[1:5]) with userdefined function

Question

I am calculating the maximum climatalogical water deficit from this study/scientific publication, using the Rcode that the authors publsihed with the study available here

It is a very simple code in which the input is monthly rainfall rasters, which I have downloaded from the Terraclimate product. After subtracting 100 mm/month (evapotranspiration) from each rainfall raster (month), a stack is created and the following function is used in 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)

However, after I run the line starting cwd I get the error Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : cannot use this function

Why do I have this error? How do I fix it?

I have looked at similar posts-1, 2,and have used na.rm=TRUE etc but I still get the same error. I also see some posts that say that the calc() is not very good, but I do not know how to solve this problem or find an alternative as I am using the code published by authors of the study I am replicating.

EDIT- I used a vector of 100 random numbers (with no NAs) to test the mcwd.f function and it works. So I think it is to do with having NA pixels, which I thought should get solved if I use na.rm=TRUE in the calc(), but it does not solve it.

Robert Hijmans · Accepted Answer · 2020-06-18T02:20:40.150

Always include a minimal, self-contained, reproducible example when you ask an R question. You say your function "works", but you do not provide input and expected output data.

I get this with your function

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

Here is a much more concise version of your function

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

And a variation that helps explain the next one

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

A vectorized version (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

Now use the function with a RasterStack and calc

Example data:

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

Use calc:

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

That works for both functions (and also for your mcwd.f). Try with NA values

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

Also works.

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

But not for your function

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

And you can see why when you do

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

To make you function work (but you should not because it is inefficient) you could add a first line like this

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

Under the assumption that when there is one NA, all values in the cell are, or should become, NA

raster::calc() error in calctest(x[1:5]) with userdefined function

1 Answers1