Connect all points with lines drawn between the closest pairs

Question

Take a set of 10 points, like this.

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"))

df.sf
# A tibble: 10 × 2
   component_number                 geometry
              <dbl> <POINT [US_survey_foot]>
 1               51       (2529694 437803.2)
 2               51       (2528863 436123.9)
 3               51       (2528991 436854.9)
 4               51       (2529139 436573.1)
 5               51       (2529133 436834.5)
 6               51       (2529134 437094.4)
 7               51       (2529134 437354.5)
 8               51       (2529193 437824.1)
 9               51       (2529456 437147.3)
10               51       (2529456 437147.3)

The points look like so: plot(df.sf)

How can I use the {{sf}} package to draw a set of each lines connecting each point to its nearest two neighbors?

Answer 1

If it doesn't have to be just sf, there are packages that can help with nearest neighbours, e.g. {spdep}:

library(tibble)
library(spdep)
library(sf)
library(ggplot2)

# remove duplicates
df.sf <- unique(df.sf)

# 2 nearest neighbours, returned knn object includes a matrix with indices:
knn <- knearneigh(df.sf, k = 2)
knn$nn
#>       [,1] [,2]
#>  [1,]    8    9
#>  [2,]    4    3
#>  [3,]    5    6
#>  [4,]    5    3
#>  [5,]    3    6
#>  [6,]    5    7
#>  [7,]    6    9
#>  [8,]    7    1
#>  [9,]    6    7

# and we can convert it to neighbours list and then to lines:
lines_sf <- nb2lines(knn2nb(knn), coords = df.sf$geometry)
lines_sf
#> Simple feature collection with 18 features and 5 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 2528863 ymin: 436123.9 xmax: 2529694 ymax: 437824.1
#> Projected CRS: NAD27 / Wisconsin South
#> First 10 features:
#>    i j i_ID j_ID wt                       geometry
#> 1  1 8    1    8  1 LINESTRING (2529694 437803....
#> 2  1 9    1    9  1 LINESTRING (2529694 437803....
#> 3  2 3    2    3  1 LINESTRING (2528863 436123....
#> 4  2 4    2    4  1 LINESTRING (2528863 436123....
#> 5  3 5    3    5  1 LINESTRING (2528991 436854....
#> 6  3 6    3    6  1 LINESTRING (2528991 436854....
#> 7  4 3    4    3  1 LINESTRING (2529139 436573....
#> 8  4 5    4    5  1 LINESTRING (2529139 436573....
#> 9  5 3    5    3  1 LINESTRING (2529133 436834....
#> 10 5 6    5    6  1 LINESTRING (2529133 436834....
plot(lines_sf$geometry)
plot(df.sf, add = TRUE)

Or work with the matrix and build our own additional geometry columns:

# build linestrings from source points(sfc) and destination index
pt_to_index_line <- function(sfc, to_idx){
    # pairwise union, returns list of multipoints
  Map(st_union, sfc, sfc[to_idx]) |> 
    # list to sfc
    st_as_sfc(crs = st_crs(df.sf)) |>
    # multipoints to linestrings
    st_cast("LINESTRING")
}

# 2 extra geometry columns for nearest neighbour links
df.sf$knn_link_1 <- pt_to_index_line(st_geometry(df.sf), knn$nn[,1])
df.sf$knn_link_2 <- pt_to_index_line(st_geometry(df.sf), knn$nn[,2])
df.sf
#> Simple feature collection with 9 features and 1 field
#> Active geometry column: geometry
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 2528863 ymin: 436123.9 xmax: 2529694 ymax: 437824.1
#> Projected CRS: NAD27 / Wisconsin South
#> # A tibble: 9 × 4
#>   component_number                 geometry                           knn_link_1
#> *            <dbl> <POINT [US_survey_foot]>        <LINESTRING [US_survey_foot]>
#> 1               51       (2529694 437803.2) (2529694 437803.2, 2529193 437824.1)
#> 2               51       (2528863 436123.9) (2528863 436123.9, 2529139 436573.1)
#> 3               51       (2528991 436854.9) (2528991 436854.9, 2529133 436834.5)
#> 4               51       (2529139 436573.1) (2529139 436573.1, 2529133 436834.5)
#> 5               51       (2529133 436834.5) (2529133 436834.5, 2528991 436854.9)
#> 6               51       (2529134 437094.4) (2529134 437094.4, 2529133 436834.5)
#> 7               51       (2529134 437354.5) (2529134 437354.5, 2529134 437094.4)
#> 8               51       (2529193 437824.1) (2529193 437824.1, 2529134 437354.5)
#> 9               51       (2529456 437147.3) (2529456 437147.3, 2529134 437094.4)
#> # ℹ 1 more variable: knn_link_2 <LINESTRING [US_survey_foot]>

ggplot(df.sf) +
  geom_sf(aes(geometry = knn_link_1, color = "knn_link_1"), linewidth = 2.5)+
  geom_sf(aes(geometry = knn_link_2, color = "knn_link_2"), linewidth = 1) +
  geom_sf() 

^{Created on 2023-11-21 with reprex v2.0.2}