在R中快速提取多边形内坐标点的高效方案咨询

阿华AIGC实验室

2026-4-21

问题背景

你提到手里有上亿级的点坐标数据集，还有数百个多边形，需要找出每个多边形包含的点。目前用sf::st_filter()处理，但数据规模太大导致速度跟不上，想找更高效的解决方案。你的可复现代码如下：

library("sf")
library("rnaturalearth")

world <- ne_countries(returnclass = "sf")
germany <- world[122, ]
plot(st_geometry(germany)) # 绘制德国边界

ex <- data.frame(lon = c(5, 12, 10),
                 lat = c(48, 50, 52))
ex_sf <- st_as_sf(ex, coords = c("lon", "lat"))
st_crs(ex_sf) <- st_crs(germany)

good_points <- st_filter(ex_sf, germany) # 筛选出德国境内的点（3个点中保留2个）

library("ggplot2")
ggplot(germany) +
  geom_sf() +
  geom_sf(data = good_points) +
  theme_bw()

高效优化方案

针对超大规模空间数据的点面匹配，我给你推荐几个经过实践验证的提速思路：

1. 强制启用空间索引（`sf`内优化）

sf的空间索引默认会自动创建，但有时候在大规模数据下可能没充分发挥作用。手动构建并启用索引，能让st_filter()先通过索引快速缩小候选点范围，再做精确的空间判断：

# 给点数据构建空间索引，同时可设置精度减少计算压力
st_geometry(ex_sf) <- st_sfc(st_geometry(ex_sf), precision = 1000)
ex_sf <- st_set_geometry(ex_sf, st_geometry(ex_sf)) # 确保索引生效

# 执行筛选时指定 predicate，进一步明确计算逻辑
good_points <- st_filter(ex_sf, germany, .predicate = st_intersects)

小提示：如果你的多边形数量也很多，同样可以给多边形数据集构建空间索引，双向优化能大幅减少计算量

2. 切换到`terra`包（更高效的空间引擎）

terra是新一代的空间数据处理包，在处理大规模矢量/栅格数据时性能远超sf，尤其是批量运算场景：

library(terra)

# 将sf格式转换为terra的vect格式
germany_vect <- vect(germany)
ex_vect <- vect(ex_sf)

# 直接用索引语法提取多边形内的点，速度极快
good_points_vect <- ex_vect[germany_vect, ]

# 如果需要转回sf格式继续处理
good_points_sf <- st_as_sf(good_points_vect)

3. 分块批量处理（解决内存瓶颈）

如果数据大到内存无法一次性加载，就把点数据拆分成若干小块，逐个处理后再合并结果：

# 根据内存情况调整分块数量，这里分成10块
chunk_num <- 10
chunks <- split(ex_sf, cut(seq_len(nrow(ex_sf)), chunk_num))

# 批量处理每个分块
good_points_list <- lapply(chunks, function(chunk) {
  st_filter(chunk, germany)
})

# 合并所有分块的结果
good_points <- do.call(rbind, good_points_list)

4. 并行计算+分块（最大化利用多核CPU）

结合并行计算框架，把分块处理的任务分配到多个CPU核心上，进一步缩短处理时间：

library(future.apply)

# 开启多会话并行（核心数可手动指定，比如workers = 4）
plan(multisession)

# 并行处理分块
good_points_list <- future_lapply(chunks, function(chunk) {
  st_filter(chunk, germany)
})

# 合并结果并关闭并行
good_points <- do.call(rbind, good_points_list)
plan(sequential)