Spatial join between two networks

[luukvdmeer]

Is your feature request related to a problem? Please describe. In tidygraph, the graph_join function can join two graphs A and B, resulting in a new graph C where nodes that were shared between A and B, only occur once in C. That is, identical nodes of the two input graphs are detected as being identical, and merged accordingly. Additionally, it updates the edges in the joining graph so they match the new indexes of the nodes in the resulting graph.

Describe the solution you'd like A function st_graph_join that does exactly what graph_join does, but based on a spatial join. That is, it will detect when nodes of A and B have the same coordinates, merge them accordingly, and update the edges.

[luukvdmeer]

I moved this to milestone 2, since it is about extending the functionalities of sf and tidygraph, beyond their current state, into the domain of spatial networks.

[luukvdmeer]

Moved this to Milestone 3 since we need some more time to find a good solution

[luukvdmeer]

Heads-up @loreabad6 @Robinlovelace @agila5 On develop there is now a function st_network_join. I am not yet sure of performance on larger networks. Lets show how it works first:

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

# First, create two networks.
n11 = st_point(c(0,0))
n12 = st_point(c(1,1))
e1 = st_sfc(st_linestring(c(n11, n12)), crs = 4326)

n21 = n12
n22 = st_point(c(0,2))
e2 = st_sfc(st_linestring(c(n21, n22)), crs = 4326)

n31 = n22
n32 = st_point(c(-1,1))
e3 = st_sfc(st_linestring(c(n31, n32)), crs = 4326)

n41 = st_point(c(0.5,0.5))
n42 = st_point(c(-1,2))
e4 = st_sfc(st_linestring(c(n41, n42)), crs = 4326)

n51 = n42
n52 = n22
e5 = st_sfc(st_linestring(c(n51, n52)), crs = 4326)

net1 = as_sfnetwork(st_as_sf(c(e1,e2,e3)))
net2 = as_sfnetwork(st_as_sf(c(e4, e5)))

par(mar = c(1,1,1,1))
plot(net1, cex = 2)
plot(net2, col = "red", add = T)

# As we can see:
# --> net1 has 4 nodes and 3 edges
# --> net2 has 3 nodes and 2 edges
# --> The networks share 1 common node
# --> A node in net2 that intersects an edge in net1 
# --> There is an edge-crossing between net1 and net2

# The basic st_network_join does the following:
# --> Full spatial join on the nodes, based on st_equals
# --> Bind edges
# --> Update from and to columns of edges accordingly

# Notice that:
# --> The 7 nodes of net1 and net2 together are joined to 6
# --> Because one node is shared spatially
# --> Edges are just binded, hence 5 resulting edges
j1 = st_network_join(net1, net2)
j1
#> # An sfnetwork with 6 nodes and 5 edges
#> #
#> # CRS:  EPSG:4326 
#> #
#> # A rooted tree with spatially explicit edges
#> #
#> # Node Data:     6 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: -1 ymin: 0 xmax: 1 ymax: 2
#>             x
#>   <POINT [°]>
#> 1       (0 0)
#> 2       (1 1)
#> 3       (0 2)
#> 4      (-1 1)
#> 5   (0.5 0.5)
#> 6      (-1 2)
#> #
#> # Edge Data:     5 x 3
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: -1 ymin: 0 xmax: 1 ymax: 2
#>    from    to                x
#>   <int> <int> <LINESTRING [°]>
#> 1     1     2       (0 0, 1 1)
#> 2     2     3       (1 1, 0 2)
#> 3     3     4      (0 2, -1 1)
#> # … with 2 more rows
plot(j1)

# We can set some additional options.

# Set blend_nodes to TRUE. This will:
# --> Split edges in x that intersect with a node in y
# --> Split edges in y that intersect with a node in x
# Notice that the plot looks the same, but:
# --> The resulting network has an extra edge.
# --> Since edge 1 in net1 is splitted by node1 of net2.
j2 = st_network_join(net1, net2, blend_nodes = T)
#> Warning: Although coordinates are longitude/latitude, st_blend assumes that they
#> are planar
#> Warning: st_blend assumes attributes are constant over geometries
j2
#> # An sfnetwork with 6 nodes and 6 edges
#> #
#> # CRS:  EPSG:4326 
#> #
#> # A directed acyclic simple graph with 1 component with spatially explicit edges
#> #
#> # Node Data:     6 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: -1 ymin: 0 xmax: 1 ymax: 2
#>             x
#>   <POINT [°]>
#> 1       (0 0)
#> 2       (1 1)
#> 3       (0 2)
#> 4      (-1 1)
#> 5   (0.5 0.5)
#> 6      (-1 2)
#> #
#> # Edge Data:     6 x 3
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: -1 ymin: 0 xmax: 1 ymax: 2
#>    from    to                x
#>   <int> <int> <LINESTRING [°]>
#> 1     5     2   (0.5 0.5, 1 1)
#> 2     1     5   (0 0, 0.5 0.5)
#> 3     2     3       (1 1, 0 2)
#> # … with 3 more rows
plot(j2)

# Set blend_crossings to TRUE. This will:
# --> Split edges in x where they cross an edge in y
# --> Split edges in y where they cross an edge in x
# Notice that:
# --> The resulting network has two extra edges.
# --> Since edge 3 in net1 is splitted by edge1 of net2
# --> And edge 1 in net2 is splitted by edge 3 of net2
j3 = st_network_join(net1, net2, blend_crossings = T)
#> Warning: Although coordinates are longitude/latitude, st_blend assumes that they
#> are planar

#> Warning: st_blend assumes attributes are constant over geometries
j3
#> # An sfnetwork with 7 nodes and 7 edges
#> #
#> # CRS:  EPSG:4326 
#> #
#> # A directed simple graph with 1 component with spatially explicit edges
#> #
#> # Node Data:     7 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: -1 ymin: 0 xmax: 1 ymax: 2
#>             x
#>   <POINT [°]>
#> 1       (0 0)
#> 2       (1 1)
#> 3       (0 2)
#> 4      (-1 1)
#> 5  (-0.5 1.5)
#> 6   (0.5 0.5)
#> # … with 1 more row
#> #
#> # Edge Data:     7 x 3
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: -1 ymin: 0 xmax: 1 ymax: 2
#>    from    to                x
#>   <int> <int> <LINESTRING [°]>
#> 1     1     2       (0 0, 1 1)
#> 2     2     3       (1 1, 0 2)
#> 3     3     5  (0 2, -0.5 1.5)
#> # … with 4 more rows
plot(j3)

# Of course we can also blend nodes AND crossings.
j4 = st_network_join(net1, net2, blend_nodes = T, blend_crossings = T)
#> Warning: Although coordinates are longitude/latitude, st_blend assumes that they
#> are planar

#> Warning: st_blend assumes attributes are constant over geometries
j4
#> # An sfnetwork with 7 nodes and 8 edges
#> #
#> # CRS:  EPSG:4326 
#> #
#> # A directed simple graph with 1 component with spatially explicit edges
#> #
#> # Node Data:     7 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: -1 ymin: 0 xmax: 1 ymax: 2
#>             x
#>   <POINT [°]>
#> 1       (0 0)
#> 2       (1 1)
#> 3       (0 2)
#> 4      (-1 1)
#> 5   (0.5 0.5)
#> 6  (-0.5 1.5)
#> # … with 1 more row
#> #
#> # Edge Data:     8 x 3
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: -1 ymin: 0 xmax: 1 ymax: 2
#>    from    to                x
#>   <int> <int> <LINESTRING [°]>
#> 1     5     2   (0.5 0.5, 1 1)
#> 2     1     5   (0 0, 0.5 0.5)
#> 3     2     3       (1 1, 0 2)
#> # … with 5 more rows
plot(j4)

# Notice that nodes are always sorted just as in x.
# This takes a little extra processing.
# If you dont care about this you can also set sort to FALSE.
j5 = st_network_join(net1, net2, blend_nodes = T, sort = F)
#> Warning: Although coordinates are longitude/latitude, st_blend assumes that they
#> are planar

#> Warning: st_blend assumes attributes are constant over geometries
j5
#> # An sfnetwork with 6 nodes and 6 edges
#> #
#> # CRS:  EPSG:4326 
#> #
#> # A directed acyclic simple graph with 1 component with spatially explicit edges
#> #
#> # Node Data:     6 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: -1 ymin: 0 xmax: 1 ymax: 2
#>             x
#>   <POINT [°]>
#> 1   (0.5 0.5)
#> 2       (1 1)
#> 3       (0 0)
#> 4       (0 2)
#> 5      (-1 1)
#> 6      (-1 2)
#> #
#> # Edge Data:     6 x 3
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: -1 ymin: 0 xmax: 1 ymax: 2
#>    from    to                x
#>   <int> <int> <LINESTRING [°]>
#> 1     1     2   (0.5 0.5, 1 1)
#> 2     3     1   (0 0, 0.5 0.5)
#> 3     2     4       (1 1, 0 2)
#> # … with 3 more rows

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

[luukvdmeer]

This is now part of the new version as st_network_join. See the join and filter vignette and the function reference.

The blend_nodes and blend_crossings arguments explained above are not part of it. They require more thought before implementing.

[kholzheu]

blend_nodes and bled_crossings are exactly what I would need for my project. They are not currently part of the code, can they be made available?