R语言编写迭代函数批量处理配对栅格计算基尼系数

你可以通过「自动匹配配对文件+封装单组处理逻辑+迭代输出」的方式实现上百组栅格的批量计算，不需要手动修改文件名重复运行代码，具体实现如下： #

1. 前置加载依赖 #

library(raster)
library(DescTools)

2. 自动匹配所有配对栅格 #

先通过文件名规则自动识别所有ntl_数字.tif和对应的pop_数字.tif，提前检查文件是否存在，避免运行中途中断：

# 替换为你本地存放栅格文件的实际文件夹路径
raster_dir <- "D:/your_raster_storage_path"

# 读取所有ntl格式栅格路径
ntl_files <- list.files(
  path = raster_dir,
  pattern = "^ntl_\\d+\\.tif$",
  full.names = TRUE
)

# 提取文件配对编号
pair_ids <- gsub(".*ntl_(\\d+)\\.tif$", "\\1", ntl_files)
# 生成对应pop栅格的路径
pop_files <- file.path(raster_dir, paste0("pop_", pair_ids, ".tif"))

# 检查缺失的配对文件
miss_pop <- !file.exists(pop_files)
if (any(miss_pop)) {
  stop(paste0("缺失以下编号对应的pop栅格：", paste(pair_ids[miss_pop], collapse = "、")))
}

3. 封装单组栅格计算函数 #

把单组处理的逻辑封装成可复用函数，处理完自动在控制台打印当前组结果：

calc_pair_gini <- function(ntl_path, pop_path, pair_id) {
  # 读取并整理ntl栅格数据
  ntl <- raster(ntl_path)
  vals_ntl <- as.data.frame(values(ntl))
  ntl_coords <- as.data.frame(xyFromCell(ntl, 1:ncell(ntl)))
  ntl_df <- cbind(ntl_coords, vals_ntl)

  # 读取pop栅格并重采样匹配ntl的范围与分辨率
  pop <- raster(pop_path)
  pop <- resample(pop, ntl, method = "bilinear")
  vals_pop <- as.data.frame(values(pop))

  # 合并、清洗数据
  block_df <- cbind(ntl_df, vals_pop)
  names(block_df)[3:4] <- c("ntl", "pop")
  block_df <- na.omit(block_df)
  block_df <- subset(block_df, select = -c(x, y))

  # 按ntl值排序后计算基尼系数
  block_df <- block_df[order(block_df$ntl), ]
  gini_val <- Gini(block_df$ntl, block_df$pop, unbiased = FALSE)

  # 控制台输出当前进度
  cat(sprintf("编号[%s]栅格组计算完成，基尼系数：%.4f\n", pair_id, gini_val))

  # 返回结构化结果
  return(data.frame(
    pair_id = pair_id,
    ntl_file = basename(ntl_path),
    pop_file = basename(pop_path),
    gini = gini_val
  ))
}

4. 批量迭代并整合结果 #

用列表存储每次循环的结果，最后一次性合并为data.frame，比逐行追加的方式效率高很多，适合上百组数据的处理：

# 初始化结果存储列表
res_list <- vector("list", length(pair_ids))

# 循环处理所有配对
for (i in seq_along(pair_ids)) {
  res_list[[i]] <- calc_pair_gini(
    ntl_path = ntl_files[i],
    pop_path = pop_files[i],
    pair_id = pair_ids[i]
  )
}

# 合并所有结果为data.frame
gini_result <- do.call(rbind, res_list)

# 可选：将结果导出为本地csv文件
write.csv(gini_result, file.path(raster_dir, "gini_result.csv"), row.names = FALSE)

可选优化：增加错误捕获 #

如果担心个别栅格损坏导致整个循环中断，可以把循环部分替换为带错误捕获的版本，自动记录处理失败的组别，不影响其他文件计算：

res_list <- vector("list", length(pair_ids))
error_log <- c()

for (i in seq_along(pair_ids)) {
  current_res <- tryCatch({
    calc_pair_gini(ntl_files[i], pop_files[i], pair_ids[i])
  }, error = function(e) {
    err_info <- sprintf("编号[%s]处理失败：%s", pair_ids[i], e$message)
    cat(err_info, "\n")
    error_log <<- c(error_log, err_info)
    return(NULL)
  })
  if (!is.null(current_res)) res_list[[i]] <- current_res
}

# 过滤失败项，合并有效结果
res_list <- res_list[!sapply(res_list, is.null)]
gini_result <- do.call(rbind, res_list)

# 输出错误日志
if (length(error_log) > 0) {
  cat("===== 处理失败记录 =====\n")
  cat(paste(error_log, collapse = "\n"))
}

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