使用R语言terra/sf工具切割GDB中重叠多边形并保留属性

一、矢量方案（优先推荐，匹配你的需求） #

你提出的「切割重叠多边形→合并属性→去重」思路完全可行，下面基于sf和terra分别给出实现代码，同时针对大数据量做优化。

1. 使用sf包实现 #

sf::st_intersection()可以自动切割所有重叠区域并保留原始属性，默认返回的结果会包含所有重叠多边形的属性列，我们只需要将同一子多边形的属性合并即可。

library(sf)
library(dplyr)

# 读取GDB数据（替换为你的文件路径和图层名）
fuel_data <- st_read("your_gdb_path.gdb", layer = "target_layer", quiet = TRUE) %>%
  st_make_valid() %>%  # 修复无效几何，避免相交报错
  st_set_precision(1e6)  # 设置空间精度，减少浮点误差导致的伪重叠

# 执行相交切割，得到所有无重叠子多边形
split_polys <- st_intersection(fuel_data)

# 合并属性：将每个子多边形对应的所有处理类型、年份去重后拼接
# 假设你的属性列名为treatment_type和treatment_year
split_polys <- split_polys %>%
  rowwise() %>%
  mutate(
    combined_treatments = paste(unique(treatment_type), collapse = "+"),
    combined_years = paste(unique(treatment_year), collapse = "+")
  ) %>%
  ungroup() %>%
  # 去重：保留几何唯一的多边形，仅保留合并后的属性和几何
  distinct(geometry, .keep_all = TRUE) %>%
  select(combined_treatments, combined_years, geometry)

2. 使用terra包实现 #

terra在处理大型矢量数据时内存效率更高，适合你的USFS大数据集：

library(terra)

# 读取数据
fuel_data <- vect("your_gdb_path.gdb", layer = "target_layer")
fuel_data <- makeValid(fuel_data)  # 修复无效几何

# 切割所有重叠区域：intersect处理自相交会生成所有子多边形
split_polys <- intersect(fuel_data, fuel_data)

# 提取原始属性，合并子多边形对应的所有处理信息
# 假设ID是原始多边形的唯一标识，treatment_type、treatment_year为目标属性
orig_treatments <- fuel_data$treatment_type
orig_years <- fuel_data$treatment_year
split_ids <- values(split_polys)[, c("ID", "ID.1")]

# 合并每个子多边形的属性
split_polys$combined_treatments <- apply(split_ids, 1, function(x) {
  paste(unique(orig_treatments[x]), collapse = "+")
})
split_polys$combined_years <- apply(split_ids, 1, function(x) {
  paste(unique(orig_years[x]), collapse = "+")
})

# 去重，保留唯一几何的多边形
unique_polys <- unique(split_polys, by = "geometry")

大数据量优化技巧 #

• 分块处理：按空间范围（比如USFS的子区域）拆分数据，处理完成后再合并结果。
• 几何简化：如果精度允许，用terra::simplify()或sf::st_simplify()减少顶点数，大幅提升处理速度。
• 并行计算：用future.apply包并行执行属性合并步骤，利用多核CPU加速。

二、栅格方案（备选） #

如果矢量处理速度无法满足需求，栅格方案可以保留所有像素的多属性信息，不会丢失重叠处理记录：

library(terra)

# 读取矢量数据
fuel_data <- vect("your_gdb_path.gdb", layer = "target_layer")
fuel_data <- makeValid(fuel_data)

# 创建参考栅格（分辨率根据业务需求调整，示例为30米）
ref_rast <- rast(fuel_data, resolution = 30)

# 生成处理类型+年份的唯一标识
fuel_data$treat_key <- paste(fuel_data$treatment_type, fuel_data$treatment_year, sep = "_")
unique_keys <- unique(fuel_data$treat_key)

# 为每个处理组合生成二进制栅格（1表示该区域有此处理，0为背景）
rast_list <- lapply(unique_keys, function(key) {
  sub_data <- fuel_data[fuel_data$treat_key == key, ]
  rasterize(sub_data, ref_rast, field = 1, background = 0)
})

# 合并为多波段栅格，每个波段对应一个处理组合
multi_band_rast <- rast(rast_list)
names(multi_band_rast) <- unique_keys

# 可选：生成属性组合栅格，每个像素存储所有处理组合的拼接结果
attribute_rast <- app(multi_band_rast, function(x) {
  active_keys <- unique_keys[x == 1]
  if (length(active_keys) == 0) return(NA)
  paste(active_keys, collapse = "+")
})

关键注意事项 #

• 必须先修复无效几何：st_make_valid()或makeValid()是避免相交、栅格化报错的前提。
• 空间精度设置：统一精度可消除浮点误差导致的伪重叠问题。
• 内存控制：处理超大数据时，优先使用terra，或采用分块策略避免内存溢出。

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