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