基于距离聚合数据框点：大样本经纬度数值求和效率优化

针对20万条经纬度点的300km范围内数值求和需求，直接用st_distance做全量距离计算会导致O(n²)的复杂度（4e10次计算），完全无法落地。以下是高效的优化方案：

核心优化思路 #

1. 转换为平面投影：WGS84（EPSG:4326）是球面坐标系，距离计算效率极低，必须转换为米单位的平面投影（如UTM分带投影），大幅提升距离计算和空间索引的效率。
2. 启用空间索引：利用空间索引将邻域查询的复杂度从O(n²)降至O(n log n)，避免不必要的距离计算。
3. 批量邻域筛选：通过st_join或st_filter结合空间索引，批量完成每个点的邻域筛选与求和。

具体实现代码 #

library(sf)
library(dplyr)

# 生成测试数据（模拟你的数据集）
set.seed(123)
test <- data.frame(
  Longitude = rnorm(2e5, 100, 10),  # 模拟东经100°左右的经度
  Latitude = rnorm(2e5, 30, 10),   # 模拟北纬30°左右的纬度
  value = runif(2e5, 0, 100)       # 待求和的数值列
)

# 转换为sf对象（WGS84坐标系）
pt_sf <- st_as_sf(test, coords = c("Longitude", "Latitude"), crs = 4326)

# 自动匹配UTM投影带（最优选择，精度高）
get_utm_crs <- function(lon, lat) {
  zone <- floor((lon + 180) / 6) + 1
  crs_code <- if (lat >= 0) paste0("326", sprintf("%02d", zone)) else paste0("327", sprintf("%02d", zone))
  as.numeric(paste0("EPSG:", crs_code))
}
utm_crs <- get_utm_crs(mean(test$Longitude), mean(test$Latitude))
pt_utm <- st_transform(pt_sf, crs = utm_crs)

# 显式创建空间索引（大幅提升邻域查询速度）
st_geometry(pt_utm) <- st_sfc(st_geometry(pt_utm), index = TRUE)

# 定义邻域半径（300km = 300000米）
radius <- 300000

# 为每个点分配唯一ID，用于分组求和
pt_utm <- pt_utm %>% mutate(id = row_number())

# 批量完成邻域筛选与求和
result <- pt_utm %>%
  st_join(
    pt_utm %>% select(value),
    join = function(x, y) st_distance(x, y) <= radius  # 自定义邻域判断逻辑
  ) %>%
  group_by(id) %>%
  summarise(sum_value = sum(value, na.rm = TRUE)) %>%
  st_drop_geometry()

# 合并结果回原始数据集
final_data <- pt_sf %>%
  mutate(id = row_number()) %>%
  left_join(result, by = "id")

进阶优化 #

如果上述方法仍有性能瓶颈，可以尝试以下方案：

• 并行计算：用furrr包开启多进程并行处理，适合多核服务器：

  library(furrr)
  plan(multisession, workers = 4)  # 根据CPU核心数调整workers
  
  sum_values <- future_map_dbl(pt_utm$id, function(i) {
    buffer <- st_buffer(pt_utm[i, ], radius)
    sum(pt_utm %>% st_filter(buffer) %>% pull(value), na.rm = TRUE)
  })
  
  pt_utm$sum_value <- sum_values
  

• 使用更轻量的包：nngeo包的st_nn函数在邻域查询上性能更优，可替代st_join；data.table结合空间索引也能进一步提升分组求和的速度。

注意事项 #

• 如果你的数据覆盖全球范围，UTM分带不适用，可改用Equal Earth投影（EPSG:6933），它是等距投影，能保证距离计算的准确性且适合全球范围。
• 绝对避免使用st_distance(pt, pt)生成全量距离矩阵，这会占用数十GB甚至上百GB内存，完全不可行。

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