R语言中如何整合栅格与数据框变量预测林分植被健康

原有流程问题说明 #

原有流程先在栅格层面计算植被健康值，再提取到林分，仅适用于所有预测变量都来自遥感栅格的场景。当存在按林分ID匹配的非栅格属性变量时，不需要强行将这类变量转换为栅格格式，调整计算顺序即可完成整合。
核心逻辑：栅格层面仅计算可从遥感影像获取的植被指数，提取到林分维度后，再通过ID匹配关联非栅格属性，最后在表格层面完成最终的植被健康值计算，避免不必要的栅格运算开销。

原有可复现代码 #

1. 依赖加载与示例数据生成 #

library(sf)
library(raster)
library(terra)
library(dplyr)

# 加载示例矢量边界
v <- vect(system.file("ex/lux.shp", package="terra"))
v <- v[c(1:12)] 
v_sf <- st_as_sf(v)

# 生成5个月份的5波段示例栅格
r <- rast(system.file("ex/elev.tif", package="terra"))
r <- rep(r, 5) * 1:5
names(r) <- paste0("band", 1:5)
ras_list <- list(r,r,r,r,r)

# 生成10个林分样点
pnts <- st_sample(v_sf, size = 10, type = "random")
pnts<- as_Spatial(pnts)

2. 原栅格层面植被健康计算逻辑 #

vis <- list()
for (i in seq_along(ras_list)) {
  b <- ras_list[[i]]
  # 仅含栅格VI的原始公式
  vis[i] <- 1.23 + 0.45*((b[[4]] + b[[3]] - b[[1]]) / (b[[4]] + b[[3]])) - 0.67*(b[[1]] * b[[3]] - b[[4]])
}

3. 原林分像元值提取逻辑 #

vi_vals <- list()
for (i in 1:length(vis)) { 
  n <- raster(vis[[i]])
  vi_vals[[i]] <- raster::extract(n, pnts, method = "bilinear") 
}

调整后的实现步骤 #

• 第一步：为林分空间对象添加唯一stand_id，确保和非栅格属性表中的ID一一对应
• 第二步：栅格层面仅计算两类植被指数，不提前计算最终健康值，减少重复运算
• 第三步：提取每个月份的VI值到林分，同时保留林分ID、添加月份标识
• 第四步：合并所有月份的提取结果，通过stand_id关联非栅格的林分属性列
• 第五步：在数据框层面代入包含所有变量（栅格VI+非栅格林分属性）的预测公式，批量计算植被健康值

修正后可运行代码 #

# --------------------------
# 1. 预处理：准备林分属性与空间ID
# --------------------------
# 模拟3列非栅格林分属性，实际使用时替换为自己的属性表即可
stand_attr <- data.frame(
  stand_id = 1:10, # 与林分空间对象一一对应
  stand_age = sample(10:50, 10, replace = T),
  soil_depth = runif(10, 30, 100),
  stand_density = runif(10, 800, 2200)
)
# 给林分空间点加上唯一ID
pnts$stand_id <- 1:10

# --------------------------
# 2. 栅格VI计算与值提取
# --------------------------
monthly_result <- list()
for (i in seq_along(ras_list)) {
  b <- ras_list[[i]]
  # 仅计算两类植被指数
  vi1 <- (b[[4]] + b[[3]] - b[[1]]) / (b[[4]] + b[[3]])
  vi2 <- b[[1]] * b[[3]] - b[[4]]
  vi_stack <- c(vi1, vi2)
  names(vi_stack) <- c("vi1", "vi2")
  
  # 提取值，直接绑定林分属性，添加月份标识
  ext_df <- terra::extract(vi_stack, vect(pnts), bind = TRUE, method = "bilinear") |>
    as.data.frame() |>
    mutate(month = i)
  
  monthly_result[[i]] <- ext_df
}
# 合并所有月份的提取结果
all_vi_df <- bind_rows(monthly_result)

# --------------------------
# 3. 关联非栅格属性，计算最终植被健康值
# --------------------------
# 按林分ID匹配属性
all_model_df <- left_join(all_vi_df, stand_attr, by = "stand_id")

# 代入包含非栅格变量的新公式计算（此处公式按需替换为自己的模型即可）
# 示例新公式：植被健康 = 1.23 + 0.45*VI1 -0.67*VI2 + 0.01*林龄 + 0.005*土层厚度 - 0.0002*林分密度
all_model_df$veg_health <- with(all_model_df,
  1.23 + 0.45*vi1 - 0.67*vi2 + 
  0.01*stand_age + 0.005*soil_depth - 0.0002*stand_density
)

注意事项 #

• 如果研究对象是面状林分而非点状样地，提取时可在terra::extract()中添加fun = mean, na.rm = T参数，先聚合得到每个林分范围内的平均VI值，后续匹配属性、计算健康值的逻辑完全一致
• 优先使用terra包的提取函数替代旧版raster::extract()，计算速度更快，bind = TRUE参数可直接保留矢量对象的属性字段，减少手动匹配ID的步骤
• 不需要将林分属性转换为栅格做栅格层面运算，该方式会生成大量冗余数据，研究区范围较大时内存占用会显著升高，表格层面的向量化运算效率远高于栅格运算

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