如何结合terra包的lapp与uniroot函数实现栅格计算？

尝试用terra的lapp函数结合uniroot对栅格像元求解方程时，出现错误：Error in f(lower, ...) : argument "I" is missing, with no default.

原错误代码 #

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)

错误原因分析 #

1. 参数顺序错误：uniroot要求第一个参数是待求解的变量（即t），但原函数把t放在中间位置，导致uniroot调用时参数错位，无法识别I等参数。
2. 变量名不一致：函数内部使用Ks（数值示例中的变量名），但栅格变量是Ksat；同时Theta_r应为小写的theta_r，大小写不匹配导致变量未定义。
3. lapp调用方式错误：lapp需要接收一个处理单个像元的函数，不能直接传递uniroot，需将uniroot包装在自定义函数中，同时处理EC和I的多组组合。
4. 固定interval不合理：数值示例中使用动态计算的max_guess作为interval上限，原代码用固定值c(0.001,100)，可能导致求解失败或结果不准确。

修正方案及完整代码 #

library(terra)

# 读取并预处理栅格
Psi_water <- rast("SUB100_120/hb_SUB100_120.tif") * 10 * -1
names(Psi_water) <- "Psi_water"
Ksat <- abs(rast("SUB100_120/ksat_SUB100_120.tif") * 240)
names(Ksat) <- "Ksat"
Lambda <- rast("SUB100_120/lambda_SUB100_120.tif")
names(Lambda) <- "Beta"
theta_r <- rast("SUB100_120/theta_r_SUB100_120.tif")
names(theta_r) <- "theta_r"
theta_s <- rast("SUB100_120/theta_s_SUB100_120.tif")
names(theta_s) <- "theta_s"

# 生物物理参数
Tp <- 6.5         
EC50 <- 8.8493    
p <- 3              
Psi_Root <- -6000  
b <- 10           

# 盐分与灌溉序列
EC_seq <- seq(0.5, 6, 1)  
I_seq <- seq(0.4, 12.2, by = 3.86)

# 修正后的目标函数：t作为第一个参数（uniroot要求的待解变量）
fct <- function(t, Psi_water, Ksat, Lambda, theta_s, theta_r, EC, I, b, Tp, EC50, p, Psi_Root) { 
    Eta <- 2 + 3 * Lambda
    Delta <- Eta / Lambda
    # 修正变量名Ks→Ksat
    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)
    # 修正变量名Ks→Ksat，Theta_r→theta_r
    Denom_inner1 <- ((I - t)/Ksat)^(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
}

# 包装函数：处理单个像元的所有EC和I组合
solve_pixel <- function(Psi_water, Ksat, Lambda, theta_s, theta_r) {
    # 计算每个I对应的max_guess（参考数值示例逻辑）
    Eta <- 2 + 3 * Lambda
    max_guess <- I_seq - Ksat * (Psi_water / Psi_Root)^Eta
    # 初始化结果向量
    res <- numeric(length(EC_seq) * length(I_seq))
    idx <- 1
    for (ec in EC_seq) {
        for (i in seq_along(I_seq)) {
            # 设置interval下限为0.001*Tp，上限为计算得到的max_guess
            interval <- c(0.001 * Tp, max_guess[i])
            # 处理可能的求解失败（比如interval无效）
            tryCatch({
                opt <- uniroot(fct, interval = interval, 
                               Psi_water = Psi_water, Ksat = Ksat, Lambda = Lambda, 
                               theta_s = theta_s, theta_r = theta_r, 
                               EC = ec, I = I_seq[i], b = b, Tp = Tp, 
                               EC50 = EC50, p = p, Psi_Root = Psi_Root)
                res[idx] <- opt$root
            }, error = function(e) {
                res[idx] <- NA  # 求解失败时赋值NA
            })
            idx <- idx + 1
        }
    }
    res
}

# 合并栅格为SpatRasterDataset
layers <- sds(Psi_water, Ksat, Lambda, theta_s, theta_r)

# 使用lapp处理每个像元，生成结果栅格
result <- lapp(layers, fun = solve_pixel)

# 为结果栅格命名（对应EC和I的组合）
result_names <- expand.grid(I = I_seq, EC = EC_seq) |> 
    apply(1, function(x) paste0("I_", x[1], "_EC_", x[2]))
names(result) <- result_names

# 查看结果
result

关键改动说明 #

• 调整fct参数顺序，将t放在第一位，符合uniroot的参数要求。
• 修正函数内部变量名不一致的问题，确保与栅格变量名匹配。
• 编写solve_pixel包装函数，循环处理所有EC和I的组合，同时计算动态的max_guess作为uniroot的interval上限，并添加异常处理避免求解失败中断程序。
• 对lapp的调用方式进行修正，传入自定义的solve_pixel函数处理每个像元。
• 为输出栅格添加清晰的命名，对应不同的EC和I组合。

内容的提问来源于stack exchange，提问作者Floyid Nicolas