How to connect edges in a network even if they dont' exactly match spatially

Question

I am using R with the following packages:

library(sf)
library(sfnetworks)
library(tibble)

I would like to create a sf network (spatial network) with some tolerance of distance between nodes

I Created these two lines :

 pts1 <-  matrix(1:4, ,2)
 ls1 <-  st_linestring(pts1) 
 pts2 <-  matrix(31:34, ,2)
 pts2[1,1] <- 1.00005 
 pts2[1,2] <- 3.00005
 ls2 <-  st_linestring(pts2)

And a network:

tibble(geometry = c(st_geometry(ls1),st_geometry(ls2))) %>% 
  st_sf() %>% 
  as_sfnetwork()

In this example I would like to connect both lines (here they are not because they don't share a common point), Is it possible to specify a tolerance regarding the distance between nodes?

Answer 1

Posting my GH anwers also here:

So point (1,3) and point (1.0005, 3.0005) should be treated as being the same node in the network right? So these two points are in the network represented by a single node. To me that raises the question: what coordinates should this combined node have? (1,3), (1.0005, 3.0005), or something in between? I would suggest to pre-process the data accordingly, before creating a network. For example:

# packages
library(sf)
#> Linking to GEOS 3.7.1, GDAL 2.2.2, PROJ 4.9.2
library(sfnetworks)
library(magrittr)

# data
pts1 <-  matrix(1:4, 2)
ls1 <-  st_linestring(pts1) 
pts2 <-  matrix(31:34, ,2)
pts2[1,1] <- 1.00005 
pts2[1,2] <- 3.00005
ls2 <-  st_linestring(pts2)
obj <- st_sf(geometry = c(st_geometry(ls1), st_geometry(ls2)))

# round coordinates
st_geometry(obj) = st_geometry(obj) %>%
  lapply(function(x) round(x, 0)) %>%
  st_sfc(crs = st_crs(obj))

obj
#> Simple feature collection with 2 features and 0 fields
#> geometry type:  LINESTRING
#> dimension:      XY
#> bbox:           xmin: 1 ymin: 3 xmax: 32 ymax: 34
#> CRS:            NA
#>                  geometry
#> 1   LINESTRING (1 3, 2 4)
#> 2 LINESTRING (1 3, 32 34)

# create network
as_sfnetwork(obj)
#> # An sfnetwork with 3 nodes and 2 edges
#> #
#> # CRS:  NA 
#> #
#> # A rooted tree with spatially explicit edges
#> #
#> # Node Data:     3 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: 1 ymin: 3 xmax: 32 ymax: 34
#>   geometry
#>    <POINT>
#> 1    (1 3)
#> 2    (2 4)
#> 3  (32 34)
#> #
#> # Edge Data:     2 x 3
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: 1 ymin: 3 xmax: 32 ymax: 34
#>    from    to     geometry
#>   <int> <int> <LINESTRING>
#> 1     1     2   (1 3, 2 4)
#> 2     1     3 (1 3, 32 34)

^{Created on 2020-08-30 by the reprex package (v0.3.0)}

Another option would be to keep the original nodes, but draw extra edges between those nodes that are within distance x from each other. A fast and dirty implementation, see below.

Note I use some internal sfnetworks functions, maybe we can export some of those. Note also that this fast implementation does not consider attributes, but only the geometry column, and will duplicate original edges when they are themselves of a distance < x. It is just a showcase, tuning will be needed to make it more useful ;) The idea, however, is the same: first pre-process the data into the desired format, and only then create the network.

# packages
library(sf)
#> Linking to GEOS 3.7.1, GDAL 2.2.2, PROJ 4.9.2
library(sfnetworks)
library(magrittr)
library(purrr)
#> 
#> Attaching package: 'purrr'
#> The following object is masked from 'package:magrittr':
#> 
#>     set_names

# data
pts1 <-  matrix(1:4, 2)
ls1 <-  st_linestring(pts1) 
pts2 <-  matrix(31:34, ,2)
pts2[1,1] <- 1.00005 
pts2[1,2] <- 3.00005
ls2 <-  st_linestring(pts2)
obj <- st_sf(geometry = c(st_geometry(ls1), st_geometry(ls2)))

connect_close_nodes = function(x, threshold) {
  # get boundary points of the edges
  nodes = sfnetworks:::get_boundary_points(obj)
  # compute distance matrix with these nodes
  dist_mat = st_distance(nodes)
  # connect those boundary points that are closer than the threshold distance
  connections = which(dist_mat < threshold, arr.ind = TRUE) %>%
    apply(1, function(x) if (x[1] != x[2]) {sfnetworks:::points_to_line(nodes[x[1],], nodes[x[2],])}) %>%
    compact() %>%
    reduce(c)
  # combine the original edges with the newly created connections
  c(st_geometry(x), connections)
}

obj = st_as_sf(connect_close_nodes(obj, threshold = 0.0001))
obj
#> Simple feature collection with 4 features and 0 fields
#> geometry type:  LINESTRING
#> dimension:      XY
#> bbox:           xmin: 1 ymin: 3 xmax: 32 ymax: 34
#> CRS:            NA
#>                                x
#> 1          LINESTRING (1 3, 2 4)
#> 2 LINESTRING (1.00005 3.00005...
#> 3 LINESTRING (1.00005 3.00005...
#> 4 LINESTRING (1 3, 1.00005 3....

as_sfnetwork(obj)
#> # An sfnetwork with 4 nodes and 4 edges
#> #
#> # CRS:  NA 
#> #
#> # A directed simple graph with 1 component with spatially explicit edges
#> #
#> # Node Data:     4 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: 1 ymin: 3 xmax: 32 ymax: 34
#>                   x
#>             <POINT>
#> 1             (1 3)
#> 2             (2 4)
#> 3 (1.00005 3.00005)
#> 4           (32 34)
#> #
#> # Edge Data:     4 x 3
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: 1 ymin: 3 xmax: 32 ymax: 34
#>    from    to                        x
#>   <int> <int>             <LINESTRING>
#> 1     1     2               (1 3, 2 4)
#> 2     3     4 (1.00005 3.00005, 32 34)
#> 3     3     1   (1.00005 3.00005, 1 3)
#> # … with 1 more row

^{Created on 2020-08-30 by the reprex package (v0.3.0)}