R中大规模数据集地理距离快速计算及邻点统计方法

现有两个包含经纬度字段的空间点位数据集：

• 大数据集（记为df1）：约2000万条观测记录
• 小数据集（记为df2）：共3.6万条观测记录
  需求为统计df1中每个点位周边200米范围内，归属df2的观测点总数。原有实现基于geosphere包的distm函数编写逐行循环，运算效率极低，运行超过24小时仍无法完成计算，原有代码如下：

library(geosphere)

#dataset1 (large)
df1 <- data.frame(longitude = c(-77.14239, -77.10750, -77.14239, -77.01797, -77.17203, -77.47230, -77.26490, -77.02824, -76.96993, -77.03185),
latitude = c(38.80575, 38.87987, 38.80575, 38.90425, 38.77076, 38.98140, 38.92800, 38.90436, 38.84569, 38.92080))

#dataset2 (small)
df2 <- data.frame(longitude = c(-75.34186, -123.59649, -108.20089, -115.16004, -87.62970),
                 latitude = c(40.11899, 44.38151, 36.71881, 36.22207, 41.71438))

# 原有逐行循环计算逻辑
for(i in 1:nrow(df1)){
  df1$n_within200m[i] <- sum(as.numeric(distm(cbind(df1$longitude[i], df1$latitude[i]),
 cbind(df2$longitude, df2$latitude)) < 200))}

• 计算量级过大：原有逻辑逐行遍历df1点位，每次都对df2全部3.6万个点计算球面距离，总计算量达到2000万*3.6万=7.2e11次，属于完全无优化的暴力计算
• 无空间剪枝：没有提前过滤经纬度差明显超过200米范围的无效点，绝大多数距离计算完全没有必要
• 循环额外开销：原有for循环未预分配结果向量内存，逐行修改数据框的操作会带来大量额外性能损耗

核心优化思路是引入空间索引，先通过经纬度范围快速筛除不可能在200米范围内的点位，仅对少量候选点计算精确距离，将总计算量降低3~4个数量级。

方案1：sf包空间索引匹配（最易实现，兼容性好） #

sf包的空间距离计算函数内置R树空间索引，不需要手动实现剪枝逻辑，代码改动量极小，2000万点位规模在普通服务器上数小时即可跑完。

library(sf)
# 将普通数据框转换为WGS84坐标系的空间对象
sf1 <- st_as_sf(df1, coords = c("longitude", "latitude"), crs = 4326)
sf2 <- st_as_sf(df2, coords = c("longitude", "latitude"), crs = 4326)

# 直接计算每个点位200米范围内的匹配点数量，dist参数单位为米
df1$n_within200m <- lengths(st_is_within_distance(sf1, sf2, dist = 200))

如果运行时内存不足，可以将df1拆分为每块100~200万行的子集，逐块计算后合并结果即可。

方案2：RANN包KD树快速检索（速度最快） #

如果对速度有更高要求，可以使用RANN包基于KD树的近邻检索算法，先快速查找每个点周边的近邻候选点，再做距离校验，速度比sf方案再快2~3倍。

library(RANN)
# 注意：nn2默认使用欧氏距离，经纬度场景下需要将搜索半径转换为度（200米约对应0.002度，低纬度区域可适当调大避免漏匹配）
res <- nn2(
  data = df2[, c("longitude", "latitude")],
  query = df1[, c("longitude", "latitude")],
  searchtype = "radius",
  radius = 0.002,
  k = nrow(df2)
)
# 统计每个点半径范围内的匹配数量
df1$n_within200m <- rowSums(res$nn.dists < 0.0018) # 用更严格的度阈值过滤，避免欧氏距离近似带来的误差

• 所有逐点全量计算距离的逻辑在千万级点位规模下都不具备可行性，必须通过空间索引做预剪枝
• 如果对距离精度要求极高，不建议用经纬度直接算欧氏距离的方案，优先选择sf包的球面距离计算接口
• 计算前提前给结果向量分配足够内存，不要在循环中逐行增长向量或修改数据框列

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