R version of https://beta.observablehq.com/@mbostock/nearest-neighbor-graph

nearest_neighbor_graph

library(sf)
library(data.table)
library(mapview)

n = 1000

x = runif(n, 0, 2)
y = runif(n, 0, 1)

pts = st_as_sf(data.frame(id_pt = 1:n, x = x, y = y), coords = c("x", "y"))
pts_union = st_union(pts)

deltri = st_triangulate(pts_union, bOnlyEdges = TRUE)

deltri = st_as_sf(
  data.frame(
    id_ln = 1:(npts(deltri)/2), 
    geometry = st_cast(deltri, "LINESTRING")
  )
)

pts_deltri_isec = as.data.table(st_join(deltri, pts))

lut = pts_deltri_isec[, (id_ln[which.min(sf:::CPL_length(geometry))]), by = "id_pt"]

shrtst = deltri[deltri$id_ln %in% lut$V1, ]

mapview(deltri, color = "#BEBEBE") +
mapview(shrtst, color = "red") + 
  mapview(pts, cex = 2, col.regions = "black")

constrain to certain distance:

library(sf)
library(data.table)
library(mapview)

n = 1000

x = runif(n, 0, 200)
y = runif(n, 0, 100)

pts = st_as_sf(data.frame(id_pt = 1:n, x = x, y = y), coords = c("x", "y"))
pts_union = st_union(pts)

deltri = st_triangulate(pts_union, bOnlyEdges = TRUE)

deltri = st_as_sf(
  data.frame(
    id_ln = 1:(npts(deltri)/2)
    , geometry = st_cast(deltri, "LINESTRING")
  )
)
deltri$lngth = st_length(deltri$geometry)
deltri = deltri[deltri$lngth <= 3, ]

mapview(deltri, color = "red", lwd = 2, legend = FALSE) +
  mapview(pts, cex = 2, col.regions = "black", legend = FALSE)