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

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

CodePudding user response：

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.