在R中能否通过循环批量读取并分析多份Shapefile？

问题背景 #

刚接触R空间数据处理，需处理近3000个按英国OS格网命名的Shapefile（总数据量数百GB，无法一次性读取）。目前能手动处理单个郡：匹配对应OS格网的Shapefile，读取合并后裁剪到郡边界，计算覆盖面积和覆盖率，但不知道如何用循环自动化，也不清楚怎么把计算结果批量加到原始county数据框中，同时想找更高效的实现方式。

现有手动处理代码：

gc()
library(sf)
library(sp)
library(ggplot2)
library(dplyr)
library(stringr)

# 创建Shapefile路径数据框
all_shps <- data_frame(filename = list.files(path = "Polygons", pattern = "\\.shp$", recursive = T))
all_shps <- all_shps  %>%
  mutate(filepath = paste0("Polygons/", filename))

# 导入基础图层
os_grid <- st_read("OSGB_Grid_1km.shp")
county <- st_read("counties.shp")

# 关联郡与OS格网ID
county_os <- st_intersection(os_grid, county)

# 手动处理第一个郡
countyx <- county$NAME[[1]]
countyname <- county %>% filter(NAME == countyx)
target <- county_os %>% filter(NAME == countyx)
targetgrid <- target$PLAN_NO

shapes <- all_shps[all_shps$filename %in% target, ] # 原代码shapes_all修正为all_shps
shapefile_list <- lapply(shapes$filepath, read_sf)
shapefile_list_all <- do.call(rbind, shapefile_list)

intersect_1 <- st_intersection(shapefile_list_all, countyname)
coverarea <- sum(st_area(intersect_1))
countyarea <- st_area(countyname)
prop_cover <- as.numeric(coverarea/countyarea)

核心问题 #

1. 能否将上述逻辑放入循环中批量处理所有郡？
2. 有没有更高效的实现方式？

解决方案 #

1. 基础循环实现 #

直接用for循环遍历每个郡，初始化结果列后批量赋值：

# 初始化结果列，提前计算郡面积避免重复计算
county <- county %>%
  mutate(cover_area = NA_real_,
         prop_cover = NA_real_,
         county_area = as.numeric(st_area(.)))

# 循环处理每个郡
for (i in seq_len(nrow(county))) {
  # 获取当前郡的信息
  current_county <- county[i, ]
  county_name <- current_county$NAME
  
  # 获取对应OS格网ID
  target_grids <- county_os %>%
    filter(NAME == county_name) %>%
    pull(PLAN_NO)
  
  # 匹配对应的Shapefile（可根据实际文件名规则调整匹配逻辑）
  target_shapes <- all_shps %>%
    filter(str_detect(filename, paste(target_grids, collapse = "|")))
  
  # 读取并裁剪每个Shapefile（先裁剪再合并，大幅减少内存占用）
  clipped_shapes <- lapply(target_shapes$filepath, function(path) {
    shp <- read_sf(path)
    st_intersection(shp, current_county)
  })
  
  # 合并裁剪后的结果
  combined_clipped <- do.call(rbind, clipped_shapes)
  
  # 计算覆盖指标
  total_cover <- as.numeric(sum(st_area(combined_clipped)))
  cover_prop <- total_cover / current_county$county_area
  
  # 赋值回原数据框
  county$cover_area[i] <- total_cover
  county$prop_cover[i] <- cover_prop
  
  # 清理临时变量，释放内存
  rm(shp, clipped_shapes, combined_clipped)
  gc()
  
  # 打印进度
  cat("完成处理：", county_name, "（", i, "/", nrow(county), "）\n", sep = "")
}

2. 更简洁高效的purrr实现 #

用purrr映射函数替代循环，结合dplyr行处理，代码更简洁：

library(purrr)

# 定义处理单个郡的函数
process_single_county <- function(county_row) {
  county_name <- county_row$NAME
  
  # 获取对应OS格网ID
  target_grids <- county_os %>%
    filter(NAME == county_name) %>%
    pull(PLAN_NO)
  
  # 匹配Shapefile
  target_shapes <- all_shps %>%
    filter(str_detect(filename, paste(target_grids, collapse = "|")))
  
  # 读取、裁剪并合并Shapefile
  combined_clipped <- target_shapes$filepath %>%
    map(read_sf) %>%
    map(~st_intersection(.x, county_row)) %>%
    list_rbind() # dplyr 1.1.0+版本可用，替代do.call(rbind)
  
  # 计算指标
  total_cover <- as.numeric(sum(st_area(combined_clipped)))
  cover_prop <- total_cover / as.numeric(st_area(county_row))
  
  # 返回带结果的郡行
  county_row %>%
    mutate(cover_area = total_cover,
           prop_cover = cover_prop)
}

# 批量处理所有郡
county_processed <- county %>%
  rowwise() %>%
  mutate(process_single_county(cur_data())) %>%
  ungroup()

3. 关键优化建议 #

• 内存优化：
  • 读取Shapefile时用st_read(..., select = c("必要列"))，只加载需要的属性列，减少内存占用。
  • 优先对单个Shapefile裁剪再合并，避免加载全部对应Shapefile后再裁剪，大幅降低内存压力。
  • 循环/映射过程中定期清理临时变量并执行gc()，强制释放内存。
• 效率提升：
  • 提前计算所有郡的面积并存储，避免重复计算。
  • 根据文件名规则优化匹配逻辑（比如文件名严格等于格网ID+后缀，用filename %in% paste0(target_grids, ".shp")更高效）。
  • 多核机器可使用furrr包实现并行处理，大幅缩短运行时间：

    library(furrr)
    plan(multisession) # 开启并行
    
    county_processed <- county %>%
      rowwise() %>%
      mutate(future_map_dfr(cur_data(), process_single_county)) %>%
      ungroup()
    

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