Connect nodes based on proximity

[luukvdmeer]

Question via email:

I would like to build a graph from provided points, with edges between the points, if the distance is <= 100 m.

A solution:

library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(sfnetworks)

# Get some example points.
pts = st_centroid(roxel[1:10, ])
#> Warning in st_centroid.sf(roxel[1:10, ]): st_centroid assumes attributes are
#> constant over geometries of x
plot(st_geometry(pts), pch = 20)

# Calculate adjacency matrix based on distance threshold.
# I used a large threshold here in the example.
adj = st_is_within_distance(pts, dist = 250, sparse = FALSE)

# Set diagonal to FALSE.
# Unless you want loop edges at every node.
diag(adj) = FALSE

# Construct edges by defining from and to node for each edge.
from = rep(c(1:nrow(pts)), apply(adj, 1, sum))
to = do.call("c", apply(adj, 1, which))
edges = data.frame(from = from, to = to)

# Construct network from scratch.
# Add linestring geometries to the edges only if you want.
net = sfnetwork(nodes = pts, edges = edges, edges_as_lines = TRUE)
#> Checking if spatial network structure is valid...
#> Spatial network structure is valid

net
#> # A sfnetwork with 10 nodes and 30 edges
#> #
#> # CRS:  EPSG:4326 
#> #
#> # A bipartite simple graph with 2 components with spatially explicit edges
#> #
#> # Node Data:     10 × 3 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: 7.532307 ymin: 51.94496 xmax: 7.543737 ymax: 51.95566
#>   name                  type                   geometry
#>   <chr>                 <fct>               <POINT [°]>
#> 1 Havixbecker Strasse   residential (7.533591 51.95566)
#> 2 Pienersallee          secondary   (7.532307 51.95375)
#> 3 Schulte-Bernd-Strasse residential (7.532794 51.95233)
#> 4 <NA>                  path        (7.539145 51.94507)
#> 5 Welsingheide          residential (7.537643 51.94656)
#> 6 <NA>                  footway     (7.543737 51.94705)
#> # … with 4 more rows
#> #
#> # Edge Data:     30 × 3
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: 7.532307 ymin: 51.94496 xmax: 7.543737 ymax: 51.95566
#>    from    to                               geometry
#>   <int> <int>                       <LINESTRING [°]>
#> 1     1     2 (7.533591 51.95566, 7.532307 51.95375)
#> 2     2     1 (7.532307 51.95375, 7.533591 51.95566)
#> 3     2     3 (7.532307 51.95375, 7.532794 51.95233)
#> # … with 27 more rows

plot(net)

^{Created on 2022-03-22 by the reprex package (v2.0.1)}

Probably there are other (better?) solutions as well. Thoughts @Robinlovelace?

On the longer term we should think of supporting more ways to construct sfnetworks from only points. This is one example, but there are many more (see #52)

[Robinlovelace]

Probably there are other (better?) solutions as well. Thoughts @Robinlovelace?

The above looks like a solid approach to me. Minimising the amount of computation done one sf objects makes sense and that looks pretty efficient, I could think of lots of applications. No guarantee that st_within_distance() is optimal but uses GEOS behind the scenes I believe and seems to work perfectly in this case. I think it could be worthy of a function, graph_from_points() or similar?

[agila5]

Hi! My proposal:

# packages
library(tidygraph)
library(sf)
library(sfnetworks)

# data
pts = st_centroid(roxel[1:10, ])
#> Warning in st_centroid.sf(roxel[1:10, ]): st_centroid assumes attributes are
#> constant over geometries of x

# Compute adj list
(adj = st_is_within_distance(pts, dist = 250))
#> Sparse geometry binary predicate list of length 10, where the predicate
#> was `is_within_distance'
#>  1: 1, 2
#>  2: 1, 2, 3
#>  3: 2, 3
#>  4: 4, 5, 7, 8, 9, 10
#>  5: 4, 5, 7, 8
#>  6: 6, 9, 10
#>  7: 4, 5, 7, 8, 9, 10
#>  8: 4, 5, 7, 8
#>  9: 4, 6, 7, 9, 10
#>  10: 4, 6, 7, 9, 10

# build edges table
edges <- data.frame(
  from = rep(seq_along(adj), times = lengths(adj)), 
  to = unlist(adj)
)
edges <- edges[with(edges, from != to), ]

# Build sfnetworks object + add explicit edges
sfn <- sfnetwork(nodes = pts, edges = edges, edges_as_lines = TRUE)
#> Checking if spatial network structure is valid...
#> Spatial network structure is valid

par(mar = rep(0, 4))
plot(sfn)

^{Created on 2022-03-23 by the reprex package (v2.0.1)}

It is more or less the same idea as @luukvdmeer's proposal, but I think it benefits from the sparse structure of sgbp objects (i.e. the size of the adj object might grow pretty fast for large n).