加速sf包st_crop处理大规模数据集的方法及多边形重构问题

处理400万个1公顷网格的空间裁剪任务，循环调用st_crop确实会因为单次操作的 overhead 导致效率极低。下面针对你的两个问题给出具体的解决思路和代码示例：

问题1：加速空间裁剪的高效方法 #

1. 利用空间索引实现批量过滤（最推荐） #

sf内置的空间索引可以大幅提升空间查询效率，避免逐个网格循环的低效操作。你可以先给目标shapefile建立空间索引，再通过批量空间操作完成裁剪：

library(sf)
library(data.table)

# 假设你有400万个网格组成的sf对象grid_sf
# 给目标shapefile（示例中的nc）建立空间索引
st_geometry(nc) <- st_geometry(nc) %>% st_set_precision(1e3) # 可选：设置精度减少计算量
st_index(nc) <- st_sfc_index(nc)

# 方式1：直接批量获取网格与多边形的交集（自动完成裁剪）
grid_with_polygons <- st_intersection(grid_sf, nc)

# 方式2：先筛选与网格范围相交的多边形，再裁剪（适合只需要多边形结果的场景）
target_extent <- st_union(grid_sf)
filtered_polygons <- st_filter(nc, target_extent)
fast_cropped_polygons <- st_crop(filtered_polygons, target_extent)

这种向量化操作会比循环快几个数量级，完全适配百万级网格的场景。

2. 先通过边界框快速过滤再裁剪 #

如果偏好data.table的操作逻辑，可以先通过多边形的边界框快速排除完全不相交的对象，再对剩余对象做裁剪：

# 提取所有多边形的边界框为data.table
nc_bbox_dt <- data.table(
  id = 1:nrow(nc),
  xmin = st_bbox(nc)[1, ],
  xmax = st_bbox(nc)[3, ],
  ymin = st_bbox(nc)[2, ],
  ymax = st_bbox(nc)[4, ]
)

# 筛选与目标范围相交的多边形ID
target_bb <- st_bbox(c(xmin= -79, xmax=-78,ymin= 34.5, ymax= 35.5))
selected_ids <- nc_bbox_dt[xmin < target_bb[3] & xmax > target_bb[1] & ymin < target_bb[4] & ymax > target_bb[2], id]

# 仅对筛选后的多边形做裁剪
fast_crop <- st_crop(nc[selected_ids, ], target_bb)

问题2：按原始顺序重构多边形 #

你之前遇到的周长翻倍问题，核心是多边形拆分后的点顺序丢失，导致重构时形成了错误的几何结构。解决关键是保留点在原多边形中的顺序信息：

步骤1：拆分多边形时记录点的顺序 #

在将多边形拆分为点的过程中，给每个点添加point_order列，标记它在原多边形中的位置：

library(mapview)

# 获取每个多边形的点数量
nobs <- npts(nc, by_feature = T)
# 提取坐标并添加ID和点顺序
NC <- data.table(
  id = rep(1:nrow(nc), nobs),
  point_order = unlist(lapply(nobs, function(x) 1:x)), # 每个多边形的点从1开始编号
  st_coordinates(nc)[, 1:2]
)

步骤2：按顺序排序后重构多边形 #

子集筛选后，先按id和point_order排序，再转换为sf对象并构建多边形：

# 子集筛选目标范围内的点
crop_NC <- NC[X >= target_bb[1] & X < target_bb[3] & Y >= target_bb[2] & Y < target_bb[4]]
# 按ID和点顺序排序，确保点的连接顺序正确
crop_NC_sorted <- crop_NC[order(id, point_order)]

# 转回sf并重构多边形
crop_NC_sf <- st_as_sf(crop_NC_sorted, coords = c("X", "Y"), crs = st_crs(nc)) %>%
  group_by(id) %>%
  summarise() %>% # 仅按ID聚合点，不需要额外计算
  st_cast("POLYGON")

现在再对比周长，结果就会和st_crop的输出一致了：

sum(st_length(fast_crop), na.rm = T) # 1307555 [m]
sum(st_length(crop_NC_sf), na.rm = T) # 与上面结果一致

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