r语言 - 如何结合unroot和terra包中的lapp/overlay函数来执行栅格计算?

Question

我正在尝试从terra包中拉普来解决具有栅格单元和单根函数的方程。我被困在如何从结合拉普和单根函数得到最终的光栅输出。我将非常感谢任何帮助

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如下:f(lower，…)出错:参数"I"缺少，没有默认值。

#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

基于您的google驱动器的栅格数据，我创建了一个自包含的可重复的示例数据(请在将来的问题中包括类似的东西):

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")

你的函数

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
}

在另一个调用uniroot的函数中使用该函数。uniroot不是矢量化的，所以我们需要做一个循环。它也没有在搜索区间内处理NA，所以我们也需要处理它。

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
}

现在用lapp调用函数。

可以修改默认参数设置。

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 

注意，当我们使用usenames=TRUE时，变量(层)名称与函数中的参数名称完全匹配。