How can I combine uniroot and lapp/overlay function from terra package to perform raster calculation?

Question

I am trying to lapp from terra package to solve an equation with raster cells and uniroot function. I am stuck on how to get the final raster output from combining lapp and uniroot functions. I would be thankful for any kind of help

library(terra)

Psi_water <- rast("SUB100_120/hb_SUB100_120.tif");
Psi_water <- Psi_water*10*-1
names(Psi_water) <- "Psi_water"
Ksat <- rast("SUB100_120/ksat_SUB100_120.tif");
Ksat <- abs(Ksat*240)
names(Ksat) <- "Ksat"
Lambda <- rast("SUB100_120/lambda_SUB100_120.tif");
names(Lambda) <- "Beta"
theta_r1 <- rast("SUB100_120/theta_r_SUB100_120.tif")
names(theta_r1) <- "theta_r"
theta_s1 <- rast("SUB100_120/theta_s_SUB100_120.tif")
names(theta_s1) <- "theta_s"

#Biophysical parameters
Tp= 6.5         
EC50= 8.8493    
p= 3              
Psi_Root=-6000  
b= 10           
Yr0= 0

#Water Salinity (EC)
EC <- seq(0.5, 6, 1)  

#Irrigation water
I <- seq(0.4 , 12.2, by =3.86)

#function 
fct <- function (Psi_water, Ksat, Lambda, theta_s, theta_r, EC, t, I, b,  Tp, EC50, p, Psi_Root) { 
    Eta <-  2 + 3 * Lambda
    Delta <- Eta / Lambda
    Num_inner1 <- ((I-t)/Ksat)^(1/Eta)
    Num_inner2 <- abs(Psi_water)/Num_inner1
    Num_inner3 <- abs(Psi_Root)-Num_inner2
    Num_inner4 <- Num_inner3*(I-t)*b
    Num <- min(Tp,Num_inner4)
    Denom_inner1 <- ((I-t)/Ks)^(1/Delta)
    Denom_inner2 <- Theta_r+(theta_s-theta_r)*Denom_inner1
    Denom_inner3 <- EC*I*Denom_inner2
    Denom_inner4 <- EC50*(I-t)*theta_s
    Denom <-  1+(Denom_inner3/Denom_inner4)^p
    Num-Denom*t #round(out,3)
}

layers <- sds(Psi_water, Ksat, Lambda, theta_s, theta_r)

t <- lapp(layers, uniroot(fct, interval = c(0001, 100)), fun=fct)

My out is the following: Error in f(lower, ...) : argument "I" is missing, with no default.

#Example with numeric values
#Biophysical parameters
Tp= 6.5         
EC50= 8.8493    
p= 3              
Psi_Root=-6000  
b= 10           
Yr0= 0

#soil parameters
Ks= 3600
Beta= 0.55
Eta= 2.7
Theta_s= 0.41
Theta_r= 0.06
Psi_Water= -200
Delta=Eta/Beta

#Water Salinity (EC)
EC <- seq(0.5, 6, 1)  

#Irrigation water
I <- seq(0.4 , 12.2, by =3.86) 

#Optimization function
fct <- function (Tp, EC50,p,Psi_Root,Ks,b,Beta,Eta,Theta_s,Theta_r,Psi_Water,Delta,EC, I,t) { 
  Num_inner1 <- ((I-t)/Ks)^(1/Eta)
  Num_inner2 <- abs(Psi_Water)/Num_inner1
  Num_inner3 <- abs(Psi_Root)-Num_inner2
  Num_inner4 <- Num_inner3*(I-t)*b
  Num <- min(Tp,Num_inner4)
  Denom_inner1 <- ((I-t)/Ks)^(1/Delta)
  Denom_inner2 <- Theta_r+(Theta_s-Theta_r)*Denom_inner1
  Denom_inner3 <- EC*I*Denom_inner2
  Denom_inner4 <- EC50*(I-t)*Theta_s
  Denom <-  1+(Denom_inner3/Denom_inner4)^p
  out <- Num-Denom*t#round(out,3)
  return(out)
}
#max guess for upper uniroot level
max_guess=(I-Ks*(Psi_Water/Psi_Root)^Eta);max_guess

t <- matrix(NA, ncol=length(EC), nrow=length(I))
rownames(t) <- I
colnames(t) <- EC

for (j in 1:length(EC)) {
  for (i in 1:length(I)) {
    max_guess[i]
    opt <- uniroot(f=fct, Tp=Tp, EC50=EC50,p=p,Psi_Root=Psi_Root,Ks=Ks,b=b,Beta=Beta,
                   Eta=Eta,Theta_s=Theta_s,Theta_r=Theta_r,Psi_Water=Psi_Water,Delta=Delta,EC=EC[j],I=I[i], interval=c(0.001*Tp, max_guess[i]))
    
    t[i, j] = opt$root 
  }
}
t
#output
> t
             0.5        1.5        2.5        3.5        4.5        5.5
0.4   0.03010785 0.03010785 0.03010785 0.03010785 0.03010785 0.03010785
4.26  3.88991889 3.88990654 3.87785872 3.72455754 3.57576228 3.43192052
8.12  6.49502364 6.38842824 6.14818024 5.85978470 5.56571210 5.27915708
11.98 6.49935915 6.48287826 6.42364157 6.30540640 6.12903770 5.90777669

Answer 1

Based on the raster data from your google drive I created a self-contained reproducible example data (please include something like that in future questions):

r <- rast(ncol=10, nrow=10)
n <- ncell(r)
set.seed(232)
Psi_Water <- setValues(r, runif(n, -8, -1))
Ks <- setValues(r, runif(n, .1, 465))
Beta <- setValues(r, runif(n, .22, .41))
theta_r <- setValues(r, runif(n, .03, .14))
theta_s <- setValues(r, runif(n, .4, .53))
x <- c(Psi_Water, Ks, Beta, theta_r, theta_s)
names(x) <- c("Psi_Water", "Ks", "Beta", "Theta_r", "Theta_s")

Your function

fct <- function(Tp, EC50, p, Theta_s, Theta_r, Psi_Water, Psi_Root, Ks, b, Beta, Eta,  Delta, EC, I, t) { 
    Num_inner1 <- ((I-t)/Ks)^(1/Eta)
    Num_inner2 <- abs(Psi_Water)/Num_inner1
    Num_inner3 <- abs(Psi_Root)-Num_inner2
    Num_inner4 <- Num_inner3*(I-t)*b
    Num <- min(Tp,Num_inner4)
    Denom_inner1 <- ((I-t)/Ks)^(1/Delta)
    Denom_inner2 <- Theta_r+(Theta_s-Theta_r)*Denom_inner1
    Denom_inner3 <- EC*I*Denom_inner2
    Denom_inner4 <- EC50*(I-t)*Theta_s
    Denom <-  1+(Denom_inner3/Denom_inner4)^p
    Num-Denom*t
}

Use that function in another function that calls uniroot. uniroot is not vectorized, so we need to do a loop. It also does not handle NA in the search interval, so we need to deal with that as well.

unifun <- function(Psi_Water, Ks, Beta, Theta_s, Theta_r, 
    Tp=6.5, EC50=8.8493, p=3, Psi_Root=-6000, Eta= 2.7,  b= 10, Yr0= 0, EC=0.5, I=0.4) {
    n <- length(Psi_Water)
    out <- rep(NA, n)
    for (i in 1:n) {
        max_guess=(I-Ks[i]*(Psi_Water[i]/Psi_Root)^Eta)
        if (is.na(max_guess)) next

        Delta=Eta/Beta[i]

        opt <- uniroot(f=fct, Theta_r=Theta_r[i], Psi_Water=Psi_Water[i], Ks=Ks[i], Theta_s=Theta_s[i], Beta=Beta[i], Tp=Tp, EC50=EC50, p=p, Psi_Root=Psi_Root, b=b, Eta=Eta, Delta=Delta,  EC=EC, I=I, interval=c(0.001*Tp, max_guess))
        out[i] <- opt$root
    }
    out
}

Now call the function with lapp. You can change default parameter settings.

lapp(x, unifun, usenames=TRUE, I=2.5, EC=8.12)
#class       : SpatRaster 
#dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source(s)   : memory
#name        :     lyr1 
#min value   : 1.203867 
#max value   : 1.704232 

Note that, as we are using usenames=TRUE, the variable (layer) names exactly match the argument names in the function.