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.
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