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.