R语言sf包如何基于多边形边界实现多边形到最近点空间连接

现有面图层文件polygon.shp，需为每个多边形匹配对应最近点，核心需求如下：

• 计算多边形边界（而非多边形质心）到最近点的距离
• 部分多边形需匹配多个最近点，匹配数量参考图层的no_matches属性字段
• 最终输出需包含字段：id、no_matches、point_id_1、distance_1、point_id_2、distance_2、point_id_3、distance_3

原测试代码运行时报错：

setwd("~/example")
library(sf)
polygon <- st_read("polygon.shp")
point <- st_read("point.shp")

nearest = try(sf::st_nearest_feature(polygon, point))
ls = sf::st_nearest_points( polygon, point[nearest], pairwise = F)

报错信息：

Error in [.data.frame(x, i) : undefined columns selected

报错原因 #

原代码存在两个核心问题：

1. st_nearest_feature仅返回每个面对应的1个最近点索引，无法支持多最近点匹配需求；当索引存在空值、长度与面图层不匹配时，直接用[取点子集就会触发列选择错误
2. st_nearest_points设置pairwise = F会返回所有面和点的组合连线，并非一一对应的最近连线，且原逻辑没有遵循面边界计算距离的规则。

可运行代码 #

注意：如果数据为经纬度地理坐标系，st_distance默认返回大圆距离（单位为米），若需要平面距离请提前将图层转换为研究区域对应的投影坐标系。

library(sf)
library(dplyr)

setwd("~/example")
# 读取数据
polygon <- st_read("polygon.shp", quiet = TRUE)
point <- st_read("point.shp", quiet = TRUE)

# 统一坐标系，避免距离计算偏差
if (st_crs(polygon) != st_crs(point)) {
  point <- st_transform(point, st_crs(polygon))
}

# 计算所有面到所有点的最短距离：sf默认计算面边界到点的最短距离，无需单独提取质心
dist_mat <- st_distance(polygon, point)

# 初始化结果表
res <- polygon |> 
  st_drop_geometry() |> 
  select(id, no_matches) |> 
  mutate(
    point_id_1 = NA_integer_,
    distance_1 = NA_real_,
    point_id_2 = NA_integer_,
    distance_2 = NA_real_,
    point_id_3 = NA_integer_,
    distance_3 = NA_real_
  )

# 逐面匹配对应数量的最近点
for (i in seq_len(nrow(polygon))) {
  # 当前面需要匹配的点数量，最多取3个（和输出字段对应）
  match_cnt <- min(polygon$no_matches[i], 3, nrow(point))
  if (match_cnt < 1) next
  
  # 对距离排序，取最近的match_cnt个点
  current_dist <- as.numeric(dist_mat[i, ])
  top_idx <- order(current_dist)[1:match_cnt]
  
  # 填充结果字段
  for (k in seq_len(match_cnt)) {
    res[i, paste0("point_id_", k)] <- point$id[top_idx[k]]
    res[i, paste0("distance_", k)] <- current_dist[top_idx[k]]
  }
}

# 导出结果
write.csv(res, "polygon_nearest_match_result.csv", row.names = FALSE, fileEncoding = "UTF-8")

逻辑说明 #

• st_distance对面要素计算到点的距离时，默认取多边形边界上离目标点最近的位置计算距离，完全符合边界距离计算要求，不需要单独提取多边形边界做二次计算
• 代码自动适配no_matches字段的匹配数量要求，当字段值小于3时仅填充对应数量的字段，剩余字段保留空值；当字段值大于3时最多返回3个最近点，和要求的输出字段结构对齐
• 预计算距离矩阵的方式比逐要素调用最近邻函数效率更高，适合万级以内要素量的计算场景。

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