R语言优化多组多边形重叠覆盖率计算效率的技术问询

针对千级数量多边形的重叠覆盖率计算，以下是几个高效的优化方向，替代低效的双重循环：

1. 利用空间索引+向量化聚合计算 #

sf内置空间索引机制，无需手动遍历每个多边形对。通过st_intersection直接计算所有相交对，再按原多边形分组聚合面积，最后计算覆盖率。

计算红色多边形被蓝色覆盖的百分比 #

library(sf)

# 假设red_sf、blue_sf是两组sf多边形对象，先添加唯一ID列
red_sf$id_red <- seq(nrow(red_sf))
blue_sf$id_blue <- seq(nrow(blue_sf))

# 计算所有相交的多边形对及相交面积
intersect_red_blue <- st_intersection(red_sf, blue_sf)
intersect_red_blue$intersect_area <- st_area(intersect_red_blue)

# 按红色多边形分组，求和相交面积并计算覆盖率
red_coverage <- aggregate(intersect_area ~ id_red, data = intersect_red_blue, sum)
red_coverage$total_area <- st_area(red_sf)[red_coverage$id_red]
red_coverage$coverage_pct <- as.numeric(red_coverage$intersect_area / red_coverage$total_area) * 100

# 合并回原红色多边形对象，无重叠的覆盖率设为0
red_sf <- merge(red_sf, red_coverage, by = "id_red", all.x = TRUE)
red_sf$coverage_pct[is.na(red_sf$coverage_pct)] <- 0

计算蓝色多边形被红色覆盖的百分比 #

# 按蓝色多边形分组聚合相交面积
blue_coverage <- aggregate(intersect_area ~ id_blue, data = intersect_red_blue, sum)
blue_coverage$total_area <- st_area(blue_sf)[blue_coverage$id_blue]
blue_coverage$coverage_pct <- as.numeric(blue_coverage$intersect_area / blue_coverage$total_area) * 100

# 合并回原蓝色多边形对象，无重叠的覆盖率设为0
blue_sf <- merge(blue_sf, blue_coverage, by = "id_blue", all.x = TRUE)
blue_sf$coverage_pct[is.na(blue_sf$coverage_pct)] <- 0

2. 提前过滤不相交的多边形对 #

通过st_filter先筛选出空间范围内可能相交的多边形，大幅减少后续计算量：

# 为每个红色多边形筛选候选蓝色多边形（仅保留相交的）
red_candidates <- st_filter(blue_sf, red_sf, .predicate = st_intersects)
# 仅对候选对计算交集，进一步降低计算负载

3. 简化几何降低计算复杂度 #

如果对多边形精度要求不高，可先简化几何对象，减少运算节点：

# 简化多边形（dTolerance单位与数据CRS一致，按需调整）
red_sf_simplified <- st_simplify(red_sf, dTolerance = 10)
blue_sf_simplified <- st_simplify(blue_sf, dTolerance = 10)

4. 并行加速大规模计算 #

针对超大规模数据，可结合furrr包实现多核并行计算：

library(furrr)
plan(multisession, workers = 4)  # 启用4核并行

# 并行计算每个红色多边形的覆盖百分比
red_sf$coverage_pct <- future_map_dbl(seq(nrow(red_sf)), function(i) {
  current_red <- red_sf[i,]
  intersects <- st_intersects(current_red, blue_sf, sparse = FALSE)[1,]
  if(any(intersects)) {
    intersect_area <- sum(st_area(st_intersection(current_red, blue_sf[intersects,])))
    as.numeric(intersect_area / st_area(current_red)) * 100
  } else {
    0
  }
})

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