0

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?

5
  • I can't find sfnetworks on CRAN. Did you get it from somewhere else or has it been removed from CRAN?
    – Spacedman
    Aug 5 '20 at 16:02
  • It is on github : github.com/luukvdmeer/sfnetworks
    – Felipe
    Aug 5 '20 at 19:54
  • Does your network not contain nodes? What type of graph are your trying to build (eg., saturated, Gabriel, SOI, kNN, stream, road, ...)? If you do not "need" to use the existing edges, for representing roads or streams, just use the nodes and reconstruct the edges. Aug 5 '20 at 21:07
  • In my real use case I have only edges(~100k lines), no nodes. here I recreated node to explain my problem). I try to build a stream network, and one of the goal is to calculate distances
    – Felipe
    Aug 5 '20 at 21:15
  • 1
    Hi @Felipe, I just opened a new issue in the github page of sfneworks to discuss this problem. I think that, at the moment, it's not possible to fix your problem but I would like to add this "tolerance" parameter as soon as possible.
    – agila
    Aug 28 '20 at 7:59
2

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.