CODEWITHSUNDEEP
rraster
tg110
tg110

Reputation: 411

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

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.

Upvotes: 0

Views: 826

Answers (2)

eeny
eeny

Reputation: 111

I had a similar issue using calc() to vectorize implementation of a custom function that used sens.slope(x) from the trend package on a RasterStack of inputs. Here is my solution, hope it helps others.

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

My input RasterStack included cells with NA's, as such I had built in NA handling in my custom function:

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

However, sens.slope(x) fails if all values in x are NA BUT ALSO if there is only one non-NA value in x. The NA handling above takes care of the former situation but not the latter. Updating my custom function to include the situation where only one of x is non-NA fixed my issue.

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

Here is how i figured it out, I pushed several dummy data sets into my function without calc() until I reproduced the error:

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

Upvotes: 0

Robert Hijmans
Robert Hijmans

Reputation: 47536

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

Upvotes: 0

Related Questions