如何使用sf包绘制连接每个点与其最近两个邻点的连线？

问题背景 #

现有一组10个点的sf格式空间数据，数据定义如下：

library(tidyverse)
library(sf)

df.sf <- structure(list(component_number = c(51, 51, 51, 51, 51, 51, 51, 
51, 51, 51), geometry = structure(list(structure(c(2529693.76455455, 
437803.242940758), class = c("XY", "POINT", "sfg")), structure(c(2528862.86355918, 
436123.858325239), class = c("XY", "POINT", "sfg")), structure(c(2528991.21479502, 
436854.889372002), class = c("XY", "POINT", "sfg")), structure(c(2529138.56071318, 
436573.087631819), class = c("XY", "POINT", "sfg")), structure(c(2529133.32326354, 
436834.480073507), class = c("XY", "POINT", "sfg")), structure(c(2529133.70746582, 
437094.447431574), class = c("XY", "POINT", "sfg")), structure(c(2529134.07395407, 
437354.456933641), class = c("XY", "POINT", "sfg")), structure(c(2529193.24413696, 
437824.056422966), class = c("XY", "POINT", "sfg")), structure(c(2529456.32924802, 
437147.290262866), class = c("XY", "POINT", "sfg")), structure(c(2529456.32924802, 
437147.290262866), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
"sfc"), precision = 0, bbox = structure(c(xmin = 2528862.86355918, 
ymin = 436123.858325239, xmax = 2529693.76455455, ymax = 437824.056422966
), class = "bbox"), crs = structure(list(input = "NAD27 / Wisconsin South", 
    wkt = "PROJCRS[\"NAD27 / Wisconsin South\",\n    BASEGEOGCRS[\"NAD27\",\n        DATUM[\"North American Datum 1927\",\n            ELLIPSOID[\"Clarke 1866\",6378206.4,294.978698213898,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4267]],\n    CONVERSION[\"Wisconsin CS27 South zone\",\n        METHOD[\"Lambert Conic Conformal (2SP)\",\n            ID[\"EPSG\",9802]],\n        PARAMETER[\"Latitude of false origin\",42,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8821]],\n        PARAMETER[\"Longitude of false origin\",-90,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8822]],\n        PARAMETER[\"Latitude of 1st standard parallel\",42.7333333333333,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8823]],\n        PARAMETER[\"Latitude of 2nd standard parallel\",44.0666666666667,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8824]],\n        PARAMETER[\"Easting at false origin\",2000000,\n            LENGTHUNIT[\"US survey foot\",0.304800609601219],\n            ID[\"EPSG\",8826]],\n        PARAMETER[\"Northing at false origin\",0,\n            LENGTHUNIT[\"US survey foot\",0.304800609601219],\n            ID[\"EPSG\",8827]]],\n    CS[Cartesian,2],\n        AXIS[\"easting (X)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"US survey foot\",0.304800609601219]],\n        AXIS[\"northing (Y)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"US survey foot\",0.304800609601219]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"United States (USA) - Wisconsin - counties of Adams; Calumet; Columbia; Crawford; Dane; Dodge; Fond Du Lac; Grant; Green; Green Lake; Iowa; Jefferson; Juneau; Kenosha; La Crosse; Lafayette; Manitowoc; Marquette; Milwaukee; Monroe; Ozaukee; Racine; Richland; Rock; Sauk; Sheboygan; Vernon; Walworth; Washington; Waukesha; Waushara; Winnebago.\"],\n        BBOX[42.48,-91.43,44.33,-86.95]],\n    ID[\"EPSG\",32054]]"), class = "crs"), n_empty = 0L)), row.names = c(NA, 
-10L), sf_column = "geometry", agr = structure(c(component_number = NA_integer_), levels = c("constant", 
"aggregate", "identity"), class = "factor"), class = c("sf", 
"tbl_df", "tbl", "data.frame"))

执行plot(df.sf)可查看点的分布，需求是：用sf包绘制连线，将每个点与其最近的两个邻点连接。

解决方案 #

步骤1：计算距离矩阵并筛选最近邻点 #

首先计算所有点之间的平面距离矩阵，为每个点筛选出距离最近的两个邻点（排除自身）：

# 计算点之间的距离矩阵
dist_matrix <- st_distance(df.sf)

# 为每个点提取最近的2个邻点（第1位是自身距离0，取第2、3位）
nearest_neighbors <- apply(dist_matrix, 1, function(x) order(x)[2:3])

# 整理为点对数据框
point_pairs <- data.frame(
  from = rep(1:nrow(df.sf), each = 2),
  to = as.vector(nearest_neighbors)
)

步骤2：生成连线的sf线段数据 #

基于点对构造LINESTRING类型的sf对象，同时去除重复线段（比如A→B和B→A属于同一条线，仅保留一次）：

# 逐个构造线段
lines <- apply(point_pairs, 1, function(pair) {
  st_sfc(
    st_linestring(c(st_coordinates(df.sf[pair[1], ]), st_coordinates(df.sf[pair[2], ]))),
    crs = st_crs(df.sf)
  )
})

# 转换为sf对象
lines_sf <- st_sf(geometry = do.call(c, lines))

# 去除重复线段：通过排序坐标字符串判断重复
lines_sf <- lines_sf %>%
  mutate(
    coord_key = map_chr(geometry, function(line) {
      coords <- st_coordinates(line)
      # 对线段的两个端点坐标排序后生成唯一标识
      sorted_coords <- coords[order(coords[,1], coords[,2]), ]
      paste(sapply(1:nrow(sorted_coords), function(i) paste(sorted_coords[i, ], collapse = ",")), collapse = ";")
    })
  ) %>%
  distinct(coord_key, .keep_all = TRUE) %>%
  select(-coord_key)

步骤3：绘制结果 #

将原始点和连线叠加绘制：

# 绘制红色点
plot(st_geometry(df.sf), pch = 16, col = "red", cex = 1.5, main = "点与最近两个邻点的连线")
# 叠加蓝色连线
plot(st_geometry(lines_sf), col = "blue", lwd = 1.5, add = TRUE)

说明 #

• 示例中存在重复点（第9、10个点），处理逻辑会自动过滤重复的连线。
• 距离计算采用数据自带坐标系（NAD27/Wisconsin South）下的平面距离，结果符合投影坐标系的度量规则。

内容的提问来源于stack exchange，提问作者John J.