基于R与terra的土壤数据插值：如何创建循环函数？

我有带坐标信息的Excel格式土壤分析数据，还有一份等高线Shapefile文件。想创建循环函数筛选目标数据并执行插值运算，优先用terra或未弃用的R包实现。运行下述代码时出现错误：

contorno<-vect(paste(talhao, safra,"polig.shp", sep="_"))

analise_solo<-read_excel(file.path(paste0(caminho_csv, paste(talhao,safra,"00_20.xlsx",sep="_"))),skip=0,col_names=TRUE, na=NS)
    
analise_solo_t <- vect(analise_solo, geom=c('Longitude', 'Latitude'), crs=crs4326, keepgeom=FALSE) 

analise_solo_dta <- analise_solo_t
  
d <- data.frame(geom(analise_solo_dta)[,c("x", "y")], as.data.frame(analise_solo_dta))

r <- rast(contorno, res=0.00002)
     
(variaveis_calcario =names(d))
(variaveis_calcario =c(variaveis_calcario[46:49],variaveis_calcario[52:53]))

gs2_ <- list()

for(i in 1:length(variaveis_calcario)){
  gs2_[i] <- gstat(formula=variaveis_calcario[i]~1, locations=~x+y, data=d, nmax=8, set=list(idp=2))
  
  mapa_ <- list() 
  for(j in 1:length(gs2_)){
    mapa_[j] <-interpolate(r, gs2_[j], debug.level=1)
    
  }
}

报错信息 #

Error in UseMethod("predict") :
没有适用于类为"list"的对象的'predict'方法
另外：警告信息：
In gs2_[i] <- gstat(formula = variaveis_calcario[i] ~ 1, locations = ~x + :
要替换的项数不是替换长度的倍数

1. 列表赋值错误：用gs2_[i] <- ...会把gstat对象包裹成子列表，导致后续插值时传入的是列表而非单个模型，这是警告和报错的核心原因，需改用gs2_[[i]] <- ...直接赋值。
2. 嵌套循环逻辑冗余：内层循环完全没必要，每次外层循环都会重置mapa_，且用未构建完成的gs2_循环，逻辑混乱，应直接在单次循环中完成模型构建与插值。
3. 公式构建严谨性：用as.formula()动态生成公式，避免字符直接传入可能的解析问题。

修正后的代码 #

library(terra)
library(gstat)
library(readxl)

# 读取边界矢量文件
contorno <- vect(paste(talhao, safra, "polig.shp", sep = "_"))

# 读取土壤分析Excel数据
analise_solo <- read_excel(file.path(caminho_csv, paste(talhao, safra, "00_20.xlsx", sep = "_")), 
                           skip = 0, col_names = TRUE, na = "NS")

# 转换为terra矢量对象（保留坐标系信息）
analise_solo_t <- vect(analise_solo, geom = c('Longitude', 'Latitude'), crs = crs4326)

# 提取坐标与属性数据合并为数据框
d <- data.frame(geom(analise_solo_t)[, c("x", "y")], as.data.frame(analise_solo_t))

# 创建匹配边界范围的栅格模板
r <- rast(contorno, res = 0.00002)

# 筛选目标插值变量
variaveis_calcario <- c(names(d)[46:49], names(d)[52:53])

# 初始化模型列表与插值结果列表
gs2_ <- vector("list", length(variaveis_calcario))
mapa_ <- vector("list", length(variaveis_calcario))

# 循环处理每个变量
for(i in seq_along(variaveis_calcario)){
  # 构建克里金插值模型
  gs2_[[i]] <- gstat(formula = as.formula(paste(variaveis_calcario[i], "~1")), 
                     locations = ~x + y, data = d, nmax = 8, set = list(idp = 2))
  # 执行插值运算
  mapa_[[i]] <- interpolate(r, gs2_[[i]])
  # 裁剪到边界范围内（可选，去除无效插值区域）
  mapa_[[i]] <- crop(mapa_[[i]], contorno)
}

# 给插值结果命名，方便后续调用
names(mapa_) <- variaveis_calcario

额外提示 #

• 确认crs4326与Shapefile的坐标系匹配，若不匹配需用project()函数统一坐标系。
• 确保na = "NS"是Excel中缺失值的正确表示，避免数据读入异常。
• 插值后裁剪边界可有效去除范围外的无效值，提升结果实用性。

内容的提问来源于stack exchange，提问作者Jean Agostini