针对网格点2000米范围的大型数据集过滤提速方案问询

我有一个包含13170180行的大型数据框，以及由1107650个点组成的网格，需要筛选每个网格点中心2000米范围内的数据点。当前使用的R代码预计耗时28天，急需优化提速。原代码如下：

library(geosphere)
library(dplyr)
library(pbapply)
library(raster)
library(sf)

# 定义网格分辨率（度）
grid_resolution <- 1 / 60  # 1弧分分辨率

# 创建空栅格网格
grid <- raster(xmn = xmin, xmx = xmax, ymn = ymin, ymx = ymax, 
               resolution = grid_resolution, crs = "+proj=longlat +datum=WGS84")

# 提取网格中心
grid_centers <- coordinates(grid)
grid_points <- data.frame(grid_centers)
colnames(grid_points) <- c("long_center", "lat_center")

# 筛选半径（2000米）
radius <- 2000

# 初始化结果列表
points_within_2000m_list <- vector("list", nrow(grid_points))

# 遍历每个网格点计算距离并筛选
points_within_2000m_list <- pblapply(1:nrow(grid_points), function(i) {
  target_lat <- grid_points[i, 2]
  target_long <- grid_points[i, 1]
  
  data$distance <- distHaversine(
    matrix(c(data$long, data$lat), ncol = 2),
    matrix(c(target_long, target_lat), ncol = 2)
  )
  
  points_within_2000m <- data %>%
    filter(distance <= radius)
  
  return(points_within_2000m)
})

核心优化思路：从O(N*M)到空间索引查询 #

原代码的核心问题是每个网格点都遍历全部1300万数据点，时间复杂度为O(110万*1300万)，完全无法高效运行。以下是针对性的提速方案：

1. 使用SF空间索引+缓冲区过滤（最有效） #

通过空间索引可以快速定位目标区域内的数据点，避免全局遍历：

library(sf)
library(dplyr)

# 转换为SF空间对象
data_sf <- st_as_sf(data, coords = c("long", "lat"), crs = "+proj=longlat +datum=WGS84")
grid_sf <- st_as_sf(grid_points, coords = c("long_center", "lat_center"), crs = "+proj=longlat +datum=WGS84")

# 转换为UTM投影（米为单位，距离计算更准确，需替换为数据所在区域的UTM编号）
utm_crs <- st_crs("EPSG:32650") # 示例：中国东部UTM投影
data_utm <- st_transform(data_sf, utm_crs)
grid_utm <- st_transform(grid_sf, utm_crs)

# 给数据建立空间索引（自动优化邻域查询速度）
st_geometry(data_utm) <- st_sfc(st_geometry(data_utm), index = TRUE)

# 生成网格点的2000米缓冲区
grid_buffer <- st_buffer(grid_utm, dist = 2000)
grid_buffer$grid_id <- seq_len(nrow(grid_buffer)) # 给每个网格点加ID

# 空间连接：一次性找出所有缓冲区内的数据点
matched_data <- st_join(data_utm, grid_buffer, join = st_within)

# 按网格ID拆分结果列表（对应原需求的list格式）
points_within_2000m_list <- split(matched_data, matched_data$grid_id)

# 可选：转换回原始经纬度坐标系
points_within_2000m_list <- lapply(points_within_2000m_list, function(x) {
  st_transform(x, st_crs(data_sf)) %>% 
    mutate(long = st_coordinates(geometry)[,1], lat = st_coordinates(geometry)[,2]) %>%
    st_drop_geometry() # 移除空间属性，回到普通数据框格式
})

2. 并行计算进一步提速 #

在空间索引的基础上，用并行处理拆分任务：

library(furrr)
plan(multisession, workers = parallel::detectCores() - 1) # 预留1核给系统

# 并行处理每个网格缓冲区
points_within_2000m_list <- future_map(seq_len(nrow(grid_buffer)), function(i) {
  target_buffer <- grid_buffer[i,]
  data_utm %>% 
    st_filter(st_within, target_buffer) %>%
    st_transform(st_crs(data_sf)) %>%
    mutate(long = st_coordinates(geometry)[,1], lat = st_coordinates(geometry)[,2]) %>%
    st_drop_geometry()
}, .progress = TRUE)

3. 基础内存与计算优化 #

• 移除原代码中data$distance的全局赋值：每次循环修改大对象会产生大量内存拷贝，严重拖慢速度。
• 用data.table替代dplyr做过滤：大数据量下data.table的子集操作速度远快于dplyr。
• 提前完成坐标转换：避免在循环内重复执行投影转换。

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