R语言处理1亿点空间连接（sf::st_join）内存溢出优化咨询

核心优化思路 #

• 取消循环中反复拼接结果表的操作，改用列表存储分片结果后一次性合并，避免大对象反复复制带来的内存和时间开销
• 提前预构建面要素的空间索引，减少每个分片匹配时的重复计算量
• 精简匹配逻辑，仅提取所需的区域ID字段，避免冗余属性的拷贝开销
• 可选择空间网格分片替代随机分组分片，减少每个分片需要匹配的面要素数量
• 可用并行框架实现多进程分片处理，充分利用CPU多核性能

优化后参考代码 #

library(data.table)
library(sf)
library(future.apply) # 可选，用于并行处理

# 投影坐标系27700使用GEOS计算性能更高，关闭s2
sf_use_s2(FALSE)

# 1. 预处理面要素，提前构建索引加速计算
Polygon <- st_read("https://ons-inspire.esriuk.com/arcgis/services/Administrative_Boundaries/Regions_December_2019_Boundaries_EN_BFE/MapServer/WFSServer?request=GetCapabilities&service=WFS") %>%
  select(rgn19cd, rgn19nm) %>% # 只保留需要的区域ID、名称字段，大幅减少内存占用
  st_make_valid() %>% # 确保面要素拓扑有效，避免匹配出错
  st_cast("MULTIPOLYGON")
# 提前预处理面要素，加速后续所有相交计算
Polygon <- st_prepare(Polygon)

# 2. 生成模拟点数据
DT <- data.table(
  X = runif(1e8, min = 300000, max = 500000),
  Y = runif(1e8, min = 200000, max = 400000),
  UID = 1:1e8,
  Group = rep(1:50, each = 2e6) # 单分片200万点，可根据自身内存大小调整分片粒度
)

# 3. 拆分数据集为列表
DT_Split <- split(DT, by = "Group")

# 可选：开启多进程并行，利用所有CPU核心
plan(multisession, workers = availableCores())

# 处理每个分片，结果统一存入列表
result_list <- future_lapply(DT_Split, function(dt) {
  # 转换为sf点对象
  pts_sf <- st_as_sf(dt, coords = c("X", "Y"), crs = 27700)
  # 空间连接匹配区域ID
  joined <- st_join(pts_sf, Polygon, join = st_intersects, left = TRUE)
  # 丢弃几何列返回纯表格
  setDT(joined)[, geometry := NULL][]
})

# 一次性合并所有分片结果
DT_Joined <- rbindlist(result_list)

额外性能提升技巧 #

• 若机器内存充足，可适当调大单个分片的大小，减少分片总数量，降低循环调度开销
• 本案例中面要素仅为英国一级行政区，数量极少，随机分片即可满足需求，无需额外做空间格网分片
• 可使用sfheaders::sf_point替代st_as_sf转换点对象，转换速度可提升3-5倍
• 若仅需要区域ID字段，可直接调用st_intersects返回的索引匹配ID，跳过sf对象属性绑定步骤，内存占用更低：

# 轻量化匹配写法，仅提取区域ID
joined_rgn_id <- Polygon$rgn19cd[st_intersects(pts_sf, Polygon) |> unlist()]
dt[, rgn19cd := joined_rgn_id]

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