考虑地球球面性的全球经纬度最近邻距离计算及R函数重投影扩展方案

原有实现的偏差来源非常明确：

• 代码中构造pt1、pt2等几何对象时，没有为裸坐标关联WGS84坐标系信息，sf包会默认将这类无CRS的几何对象按平面坐标规则处理
• 平面规则下计算的最近邻连线长度，本质是经纬度数值的欧氏差值，输出单位是角度，完全没有实际物理距离意义，既没有适配地球曲面特性，也达不到公里级的计算精度要求

适配逻辑说明 #

sf包原生支持椭球测地线计算能力，不需要手动做平面投影转换，只要为几何对象正确绑定WGS84坐标系，开启球面计算开关，就能直接输出符合物理意义的实际距离，从根源上规避平面投影的变形误差，尤其适配全球范围的点位距离计算场景。

修正后可运行代码 #

library(raster)
library(purrr)
library(sf)
library(dplyr)

# 构造模型输出的存在点数据
r1 <- raster(nrow=360, ncol=720)
crs(r1) <- "+proj=longlat +datum=WGS84 +no_defs"
values(r1) <- rbinom(ncell(r1), 2, 0.01)
r1_points <- rasterToPoints(r1)
r1_df <- data.frame(r1_points)
r1_presence <- r1_df %>% dplyr::filter(layer==1)
# 显式将存在点转为绑定WGS84坐标系(EPSG:4326)的sf矢量对象
presence_sf <- st_as_sf(r1_presence, coords = c("x", "y"), crs = "EPSG:4326")

# 构造调查点位数据
survey_coords <- cbind(rnorm(50) * 5 + 10, rnorm(50) + 50)
# 同样绑定WGS84坐标系，避免sf误判为平面坐标
survey_sf <- st_as_sf(data.frame(survey_coords), coords = c("X1", "X2"), crs = "EPSG:4326")

# 扩展后的距离计算函数
get_distances <- function(point_idx, target_sf, source_sf, use_geodesic = TRUE) {
  # 取出对应索引的单源存在点
  single_pt <- source_sf[point_idx, ]
  # 计算与目标调查点的最近邻连线，geodesic参数控制是否启用球面计算
  nearest_line <- st_nearest_points(single_pt, target_sf, geodesic = use_geodesic)
  # 直接输出以公里为单位的距离数值，不需要手动做单位转换
  return(as.numeric(st_length(nearest_line, unit = "km")))
}

# 批量计算所有存在点的最近邻距离
output <- map_dbl(
  .x = 1:nrow(presence_sf),
  .f = get_distances,
  target_sf = survey_sf,
  source_sf = presence_sf,
  use_geodesic = TRUE
)

函数关键配置说明 #

• 新增use_geodesic开关参数：设为TRUE时，sf将基于WGS84参考椭球的测地线规则计算最短路径，完全适配地球曲面特性；设为FALSE时可回退到传统平面计算模式
• 强制通过st_as_sf为坐标绑定EPSG:4326（WGS84经纬度）坐标系，从根源上避免sf对坐标类型的误判
• 距离输出直接指定km单位，依托sf内置的椭球参数完成精度校准，不需要手动做弧度到公里的换算

• 局部等距投影适配
  针对小范围研究区域，可以选择以研究区几何中心为原点的方位等距投影，将所有点位投影到该平面坐标系下后计算欧氏距离。这种方式小区域精度足够，但在全球尺度下会存在明显的边缘投影变形，不适合跨区域的全球点位计算。
• 专用球面距离公式直接计算
  可以采用Haversine公式、Vincenty公式直接推导两点间的球面/椭球距离，不需要依赖sf的几何计算流程。这类方案计算速度更快，适合超大规模点位的批量匹配；其中Vincenty公式针对WGS84椭球的计算精度可达毫米级，比简化的球面Haversine公式更贴近实际物理距离。
• 空间索引优化计算效率
  当点位规模较大时，可以先用st_nearest_feature提前完成所有存在点与调查点的最近邻索引匹配，再批量计算匹配点对的距离，替代逐点构造几何对象的map循环，能大幅降低计算耗时，提升批量处理效率。

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