R中按区域分组计算研究位点与采集位点的空间米制距离

安全阿岳

2026-10-4

实现方案

优先推荐使用data.table实现，该方案针对大体积数据集做了优化，按引用修改数据无需复制全量数据，内存占用低、运行速度远高于普通循环和foreach方案。

依赖安装

install.packages(c("data.table", "geosphere"))

data.table：高性能结构化数据处理包，适配千万行级数据集运算
geosphere：用于经纬度坐标转米制距离，若你的x/y已经是投影后的平面米坐标，可不需要该包

完整代码

# 加载包
library(data.table)
library(geosphere)

# 导入你的数据集，示例数据构造如下
exdb <- data.frame(zone = c(1,1,1,2,2,2),
                   site = c("study", "collect", "collect", "study", "collect", "collect"),
                   x = c(53.307726, 53.310660, 53.307089, 53.313831, 53.319087, 53.318792),
                   y = c(-6.222291, -6.217151, -6.215080, -6.214152, -6.218723, -6.215815))

# 转换为data.table格式
setDT(exdb)

# 核心分组计算逻辑
exdb[, dist := {
  # 提取当前zone唯一study位点的坐标
  std_x <- x[site == "study"][1]
  std_y <- y[site == "study"][1]
  # 计算距离，study位点直接设为0
  ifelse(site == "study", 
         0, 
         # 经纬度坐标用这个计算米制距离，输入顺序为(经度, 纬度)
         distHaversine(cbind(y, x), cbind(std_y, std_x))
         # 如果是平面米坐标，替换上面一行代码为欧氏距离计算即可：sqrt((x - std_x)^2 + (y - std_y)^2)
  )
}, by = zone]

# 输出结果查看
print(exdb)