library(sfnetworks)
library(sf)
library(tidygraph)
library(ggraph)
library(dplyr)Spatial joins and filters
The integration with {sf} and addition of several spatial network specific functions in {sfnetworks} allow to easily filter information from a network based on spatial relationships, and to join new information into a network based on spatial relationships. This vignette presents several ways to do that.
Both spatial filters and spatial joins use spatial predicate functions to examine spatial relationships. Spatial predicates are mathematically defined binary spatial relations between two simple feature geometries. Often used examples include the predicate equals (geometry x is equal to geometry y) and the predicate intersects (geometry x has at least one point in common with geometry y). For an overview of all available spatial predicate functions in {sf} and links to detailed explanations of the underlying algorithms, see here.
Spatial filters
Using st_filter
Information can be filtered from a network by using spatial predicate functions inside the sf function sf::st_filter(), which works as follows: the function is applied to a set of geometries A with respect to another set of geometries B, and removes features from A based on their spatial relation with the features in B. A practical example: when using the predicate intersects, all geometries in A that do not intersect with any geometry in B are removed.
When applying sf::st_filter() to a sfnetwork, it is internally applied to the active element of that network. For example: filtering information from a network A with activated nodes, using a set of polygons B and the predicate intersects, will remove those nodes that do not intersect with any of the polygons in B from the network. When edges are active, it will remove the edges that do not intersect with any of the polygons in B from the network.
Although the filter is applied only to the active element of the network, it may also affect the other element. When nodes are removed, their incident edges are removed as well. However, when edges are removed, the nodes at their endpoints remain, even if they don’t have any other incident edges. This behavior is inherited from {tidygraph} and understandable from a graph theory point of view: by definition nodes can exist peacefully in isolation, while edges can never exist without nodes at their endpoints. The isolated nodes that remain after filtering the edges can be easily removed using a combination of a regular dplyr::filter() verb and the tidygraph::node_is_isolated() query function.
net = as_sfnetwork(mozart, "gabriel")
ply = st_buffer(st_centroid(st_combine(mozart)), 300)
filtered_by_nodes = net |>
st_filter(ply, .pred = st_intersects)
filtered_by_edges_a = net |>
activate(edges) |>
st_filter(ply, .pred = st_intersects)
filtered_by_edges_b = net |>
activate(edges) |>
st_filter(ply, .pred = st_intersects) |>
activate(nodes) |>
filter(!node_is_isolated())ggraph(net, "sf") +
geom_edge_sf() +
geom_node_sf(size = 4) +
geom_sf(data = ply, color = "orange", fill = NA, linewidth = 1) +
theme_void()
ggraph(filtered_by_nodes, "sf") +
geom_edge_sf() +
geom_node_sf(size = 4) +
geom_sf(data = ply, color = "orange", fill = NA, linewidth = 1) +
theme_void()
ggraph(filtered_by_edges_a, "sf") +
geom_edge_sf() +
geom_node_sf(size = 4) +
geom_sf(data = ply, color = "orange", fill = NA, linewidth = 1) +
theme_void()
ggraph(filtered_by_edges_b, "sf") +
geom_edge_sf() +
geom_node_sf(size = 4) +
geom_sf(data = ply, color = "orange", fill = NA, linewidth = 1) +
theme_void()
For non-spatial filters applied to attribute columns, simply use dplyr::filter() instead of sf::st_filter().
Using spatial node and edge query functions
In {tidygraph}, filtering information from networks is done by using specific node or edge query functions inside the dplyr::filter() verb. An example was already shown above, where isolated nodes were removed from the network.
In {sfnetworks}, several spatial predicates are implemented as node and edge query functions such that you can also do spatial filtering in tidygraph style. See here for a list of all implemented spatial node query functions, and here for the spatial edge query functions. Using them makes spatial filter operations fit better into the tidy workflows of {tidygraph}. For example, we could filter edges that do not cross any other edge. The tidygraph::.E() function used in the example makes it possible to directly access the complete edges table inside verbs. Similarly, we can use tidygraph::.N() to access the nodes table and tidygraph::.G() to access the network object as a whole.
complete_net = as_sfnetwork(mozart, "complete")
filtered = complete_net |>
activate(edges) |>
filter(!edge_crosses(.E())) |>
activate(nodes) |>
filter(!node_is_isolated())ggraph(complete_net, "sf") +
geom_edge_sf() +
geom_node_sf(size = 4) +
theme_void()
ggraph(filtered, "sf") +
geom_edge_sf() +
geom_node_sf(size = 4) +
theme_void()
Besides predicate query functions, you can also use the coordinate query functions for spatial filters on the nodes. For example:
x = 4549358
filtered_by_coords = net |>
filter(node_X() > x)ggraph(net, "sf") +
geom_edge_sf() +
geom_node_sf(size = 4) +
geom_vline(xintercept = x, linewidth = 1, color = "orange") +
theme_void()
ggraph(filtered_by_coords, "sf") +
geom_edge_sf() +
geom_node_sf(size = 4) +
geom_vline(xintercept = x, linewidth = 1, color = "orange") +
theme_void()
Clipping
Filtering returns a subset of the original geometries, but leaves those geometries themselves unchanged. This is different from clipping, in which they get cut at the border of a provided clip feature. There are three ways in which you can do this: sf::st_intersection() keeps only those parts of the original geometries that lie within the clip feature, sf::st_difference() keeps only those parts of the original geometries that lie outside the clip feature, and sf::st_crop() keeps only those parts of the original geometries that lie within the bounding box of the clip feature.
Note that in the case of the nodes, clipping is not different from filtering, since point geometries cannot fall party inside and partly outside another feature. However, in the case of the edges, clipping will cut the linestring geometries of the edges at the border of the clip feature (or in the case of cropping, the bounding box of that feature). To preserve a valid spatial network structure, {sfnetworks} adds new nodes at these cut locations.
clipped = net |>
activate(edges) |>
st_intersection(ply) |>
activate(nodes) |>
filter(!node_is_isolated())
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometriesggraph(net, "sf") +
geom_edge_sf() +
geom_node_sf(size = 4) +
geom_sf(data = ply, color = "orange", fill = NA, linewidth = 1) +
theme_void()
ggraph(filtered_by_edges_b, "sf") +
geom_edge_sf() +
geom_node_sf(size = 4) +
geom_sf(data = ply, color = "orange", fill = NA, linewidth = 1) +
theme_void()
ggraph(clipped, "sf") +
geom_edge_sf() +
geom_node_sf(size = 4) +
geom_sf(data = ply, color = "orange", fill = NA, linewidth = 1) +
theme_void()
Spatial joins
Using st_join
Information can be spatially joined into a network by using spatial predicate functions inside the sf function sf::st_join(), which works as follows: the function is applied to a set of geometries A with respect to another set of geometries B, and attaches feature attributes from features in B to features in A based on their spatial relation. A practical example: when using the predicate intersects, feature attributes from feature y in B are attached to feature x in A whenever x intersects with y.
When applying sf::st_join() to a sfnetwork object, it is internally applied to the active element of that network. For example: joining information into network A with activated nodes, from a set of polygons B and using the predicate intersects, will attach attributes from a polygon in B to those nodes that intersect with that specific polygon. When edges are active, it will attach the same information but to the intersecting edges instead.
Lets show this with an example in which we first create imaginary postal code areas for the Mozart dataset.
codes = net |>
st_make_grid(n = c(2, 2)) |>
st_as_sf() |>
mutate(code = as.character(seq(1000, 1000 + n() * 10 - 10, 10)))
joined = net |>
st_join(codes, join = st_intersects)
joined#> # A sfnetwork: 17 nodes and 50 edges
#> #
#> # A bipartite simple graph with 1 component and spatially explicit edges
#> #
#> # Dimension: XY
#> # Bounding box: xmin: 4548664 ymin: 2747309 xmax: 4549589 ymax: 2748537
#> # Projected CRS: ETRS89-extended / LAEA Europe
#> #
#> # Node data: 17 × 5 (active)
#> name type website geometry code
#> <chr> <chr> <chr> <POINT [m]> <chr>
#> 1 Mozartkino cinema https://www.mo… (4549504 2747309) 1010
#> 2 Haus für Mozart theatre NA (4549003 2747376) 1000
#> 3 Mozartsteg/Rudolfskai bus_stop NA (4549589 2747507) 1010
#> 4 Mozart Denkmal artwork NA (4549387 2747514) 1010
#> 5 Mozartsteg/Rudolfskai bus_stop NA (4549491 2747551) 1010
#> 6 Mozartsteg bridge NA (4549473 2747624) 1010
#> # ℹ 11 more rows
#> #
#> # Edge data: 50 × 3
#> from to geometry
#> <int> <int> <LINESTRING [m]>
#> 1 1 3 (4549504 2747309, 4549589 2747507)
#> 2 1 4 (4549504 2747309, 4549387 2747514)
#> 3 2 4 (4549003 2747376, 4549387 2747514)
#> # ℹ 47 more rows
ggraph(net, "sf") +
geom_sf(data = codes, aes(fill = code)) +
geom_edge_sf() +
geom_node_sf(size = 4) +
theme_void()
ggraph(joined, "sf") +
geom_edge_sf() +
geom_node_sf(aes(color = code), size = 4) +
theme_void()
In the example above, the polygons are spatially distinct. Hence, each node can only intersect with a single polygon. But what would happen if we do a join with polygons that overlap? The attributes from which polygon will then be attached to a node that intersects with multiple polygons at once? In {sf} this issue is solved by duplicating such a point as much times as the number of polygons it intersects with, and attaching attributes of each intersecting polygon to one of these duplicates. This approach does not fit the network case, however. An edge can only have a single node at each of its endpoints, and thus, the duplicated nodes will be isolated and redundant in the network structure. Therefore, {sfnetworks} will only join the information from the first match whenever there are multiple matches for a single node. A warning is given in that case such that you are aware of the fact that not all information was joined into the network.
Only when you set ignore_multiple = FALSE, multiple matches will result in duplicated nodes, but these duplicates are isolated (i.e. not connected to the rest of the network). You can then use the morpher to_spatial_unique() to merge spatially duplicated nodes into one, specifying how their attributes should be combined. See here for an example.
Note that in the case of joining on the edges, multiple matches per edge are not a problem for the network structure. It will simply duplicate the edge (i.e. creating a set of parallel edges) whenever this occurs.
box = st_as_sfc(st_bbox(mozart))
two_equal_polys = st_as_sf(c(box, box)) |>
mutate(foo = c("a", "b"))
# Join on nodes gives a warning that only the first match per node is joined.
# The number of nodes in the resulting network remains the same.
net |>
st_join(two_equal_polys, join = st_intersects)
#> Warning: [1m[22m`st_join()` did not join all features.
#> [33m![39m Multiple matches were detected for some nodes, of which all but the first one
#> are ignored.
#> [36mℹ[39m If you want to add multiple matches as isolated nodes instead, set
#> `ignore_multiple` to `FALSE`.#> # A sfnetwork: 17 nodes and 50 edges
#> #
#> # A bipartite simple graph with 1 component and spatially explicit edges
#> #
#> # Dimension: XY
#> # Bounding box: xmin: 4548664 ymin: 2747309 xmax: 4549589 ymax: 2748537
#> # Projected CRS: ETRS89-extended / LAEA Europe
#> #
#> # Node data: 17 × 5 (active)
#> name type website geometry foo
#> <chr> <chr> <chr> <POINT [m]> <chr>
#> 1 Mozartkino cinema https://www.mo… (4549504 2747309) a
#> 2 Haus für Mozart theatre NA (4549003 2747376) a
#> 3 Mozartsteg/Rudolfskai bus_stop NA (4549589 2747507) a
#> 4 Mozart Denkmal artwork NA (4549387 2747514) a
#> 5 Mozartsteg/Rudolfskai bus_stop NA (4549491 2747551) a
#> 6 Mozartsteg bridge NA (4549473 2747624) a
#> # ℹ 11 more rows
#> #
#> # Edge data: 50 × 3
#> from to geometry
#> <int> <int> <LINESTRING [m]>
#> 1 1 3 (4549504 2747309, 4549589 2747507)
#> 2 1 4 (4549504 2747309, 4549387 2747514)
#> 3 2 4 (4549003 2747376, 4549387 2747514)
#> # ℹ 47 more rows
# With these settings multiple matches result in duplicated nodes.
# In this example it means we have twice the number of nodes than before.
# The duplicated nodes are isolated, i.e. not connected to any other node.
net |>
st_join(two_equal_polys, join = st_intersects, ignore_multiple = FALSE)
#> Warning: [1m[22m`st_join()` created isolated nodes.
#> [33m![39m Multiple matches were detected for some nodes, of which all but the first one
#> are added as isolated nodes to the network.
#> [36mℹ[39m If you want to ignore multiple matches instead, set `ignore_multiple` to
#> `TRUE`.#> # A sfnetwork: 34 nodes and 50 edges
#> #
#> # A bipartite simple graph with 18 components and spatially explicit edges
#> #
#> # Dimension: XY
#> # Bounding box: xmin: 4548664 ymin: 2747309 xmax: 4549589 ymax: 2748537
#> # Projected CRS: ETRS89-extended / LAEA Europe
#> #
#> # Node data: 34 × 5 (active)
#> name type website foo geometry
#> <chr> <chr> <chr> <chr> <POINT [m]>
#> 1 Mozartkino cinema https://www.mo… a (4549504 2747309)
#> 2 Haus für Mozart theatre NA a (4549003 2747376)
#> 3 Mozartsteg/Rudolfskai bus_stop NA a (4549589 2747507)
#> 4 Mozart Denkmal artwork NA a (4549387 2747514)
#> 5 Mozartsteg/Rudolfskai bus_stop NA a (4549491 2747551)
#> 6 Mozartsteg bridge NA a (4549473 2747624)
#> # ℹ 28 more rows
#> #
#> # Edge data: 50 × 3
#> from to geometry
#> <int> <int> <LINESTRING [m]>
#> 1 1 3 (4549504 2747309, 4549589 2747507)
#> 2 1 4 (4549504 2747309, 4549387 2747514)
#> 3 2 4 (4549003 2747376, 4549387 2747514)
#> # ℹ 47 more rows
# Join on edges duplicates edges that have multiple matches.
# The number of edges in the resulting network is higher than in the original.
net |>
activate(edges) |>
st_join(two_equal_polys, join = st_intersects)#> # A sfnetwork: 17 nodes and 100 edges
#> #
#> # A bipartite multigraph with 1 component and spatially explicit edges
#> #
#> # Dimension: XY
#> # Bounding box: xmin: 4548664 ymin: 2747309 xmax: 4549589 ymax: 2748537
#> # Projected CRS: ETRS89-extended / LAEA Europe
#> #
#> # Edge data: 100 × 4 (active)
#> from to geometry foo
#> <int> <int> <LINESTRING [m]> <chr>
#> 1 1 3 (4549504 2747309, 4549589 2747507) a
#> 2 1 3 (4549504 2747309, 4549589 2747507) b
#> 3 1 4 (4549504 2747309, 4549387 2747514) a
#> 4 1 4 (4549504 2747309, 4549387 2747514) b
#> 5 2 4 (4549003 2747376, 4549387 2747514) a
#> 6 2 4 (4549003 2747376, 4549387 2747514) b
#> # ℹ 94 more rows
#> #
#> # Node data: 17 × 4
#> name type website geometry
#> <chr> <chr> <chr> <POINT [m]>
#> 1 Mozartkino cinema https://www.mozartki… (4549504 2747309)
#> 2 Haus für Mozart theatre NA (4549003 2747376)
#> 3 Mozartsteg/Rudolfskai bus_stop NA (4549589 2747507)
#> # ℹ 14 more rows
For non-spatial joins based on attribute columns, simply use a join function from dplyr (e.g. dplyr::left_join() or dplyr::inner_join()) instead of sf::st_join().
Join points to their nearest node
Another network specific use-case of spatial joins would be to join information from external points of interest (POIs) into the nodes of the network. However, to do so, such points need to have exactly equal coordinates to one of the nodes. Often this will not be the case. To solve such situations, you will first need to update the coordinates of the POIs to match those of their nearest node. This can be done using st_project_on_network().
# Create a network.
n1 = st_point(c(0, 0))
n2 = st_point(c(1, 0))
net = st_sf(geometry = st_sfc(st_linestring(c(n1, n2)))) |>
as_sfnetwork()
# Create a set of POIs.
p1 = st_point(c(0, 0.2))
p2 = st_point(c(0.6, 0.2))
pois = st_sf(
poi_type = c("bakery", "butcher"),
geometry = st_sfc(p1, p2)
)
# Update coordinates of POIs to match their nearest node.
ppois = st_project_on_network(pois, net, on = "nodes")ggraph(net, "sf") +
geom_sf(
data = pois, aes(color = poi_type),
pch = 8, size = 4, show.legend = FALSE
) +
geom_edge_sf() +
geom_node_sf(size = 4) +
theme_void()
ggraph(net, "sf") +
geom_sf(
data = st_nearest_points(pois, st_combine(ppois)),
color = "grey",
linetype = 2
) +
geom_sf(
data = pois, aes(color = poi_type),
pch = 8, size = 4, show.legend = FALSE
) +
geom_edge_sf() +
geom_node_sf(size = 4) +
geom_sf(
data = ppois, aes(color = poi_type),
pch = 8, size = 4, show.legend = FALSE
) +
theme_void()
After snapping the POIs, we can use sf::st_join() as expected. Do remember that if multiple POIs are snapped to the same node, only the information of the first one is joined into the network, unless you set ignore_multiple = FALSE.
st_join(net, ppois)#> # A sfnetwork: 2 nodes and 1 edges
#> #
#> # A rooted tree with spatially explicit edges
#> #
#> # Dimension: XY
#> # Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 0
#> # CRS: NA
#> #
#> # Node data: 2 × 2 (active)
#> geometry poi_type
#> <POINT> <chr>
#> 1 (0 0) bakery
#> 2 (1 0) butcher
#> #
#> # Edge data: 1 × 3
#> from to geometry
#> <int> <int> <LINESTRING>
#> 1 1 2 (0 0, 1 0)
Blending points into a network
In the example above, it makes sense to include the information from the first POI in an already existing node. For the second POI, however, its nearest node is quite far away relative to the nearest location on its nearest edge. In that case, you might want to split the edge at that location, and add a new node to the network. For this combination process we use the metaphor of throwing the network and POIs together in a blender, and mix them smoothly together.
The function st_network_blend() does exactly that. For each POI, it finds the nearest location p on the nearest edge e. If p is an already existing node (i.e. p is an endpoint of e), it joins the information from the POI into that node. If p is not an already existing node, it subdivides e at p, adds p as a new node to the network, and joins the information from the POI into that new node. For this process, it does not matter if p is an interior point in the linestring geometry of e.
blend = st_network_blend(net, pois)
blend#> # A sfnetwork: 3 nodes and 2 edges
#> #
#> # A rooted tree with spatially explicit edges
#> #
#> # Dimension: XY
#> # Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 0
#> # CRS: NA
#> #
#> # Node data: 3 × 2 (active)
#> geometry poi_type
#> <POINT> <chr>
#> 1 (0 0) bakery
#> 2 (1 0) NA
#> 3 (0.6 0) butcher
#> #
#> # Edge data: 2 × 3
#> from to geometry
#> <int> <int> <LINESTRING>
#> 1 1 3 (0 0, 0.6 0)
#> 2 3 2 (0.6 0, 1 0)
ggraph(net, "sf") +
geom_sf(
data = pois, aes(color = poi_type),
pch = 8, size = 4, show.legend = FALSE
) +
geom_edge_sf() +
geom_node_sf(size = 4) +
theme_void()
ggraph(blend, "sf") +
geom_sf(
data = st_nearest_points(pois, st_combine(st_geometry(blend, "nodes"))),
color = "grey",
linetype = 2
) +
geom_sf(
data = pois, aes(color = poi_type),
pch = 8, size = 4, show.legend = FALSE
) +
geom_edge_sf() +
geom_node_sf(size = 4) +
theme_void()
The st_network_blend() function has a tolerance parameter, which defines the maximum distance a POI can be from the network in order to be blended in. Hence, only the POIs that are at least as close to the network as the tolerance distance will be blended, and all others will be ignored. The tolerance can be specified as a non-negative number. By default it is assumed its units are meters, but this behavior can be changed by manually setting its units with units::units().
# Update the POIs.
p3 = st_point(c(0.4, 0.3))
pois = st_sf(
poi_type = c("bakery", "butcher", "bar"),
geometry = st_sfc(p1, p2, p3)
)
blend = st_network_blend(net, pois)
blend_with_tolerance = st_network_blend(net, pois, tolerance = 0.2)ggraph(blend, "sf") +
geom_sf(
data = st_nearest_points(pois, st_combine(st_geometry(blend, "nodes"))),
color = "grey",
linetype = 2
) +
geom_sf(
data = pois, aes(color = poi_type),
pch = 8, size = 4, show.legend = FALSE
) +
geom_edge_sf() +
geom_node_sf(size = 4) +
theme_void()
ggraph(blend_with_tolerance, "sf") +
geom_sf(
data = st_nearest_points(pois, st_combine(st_geometry(blend, "nodes"))),
color = "grey",
linetype = 2
) +
geom_sf(
data = pois, aes(color = poi_type),
pch = 8, size = 4, show.legend = FALSE
) +
geom_edge_sf() +
geom_node_sf(size = 4) +
theme_void()
There are a few important details to be aware of when using st_network_blend(). Firstly: when multiple POIs have the same nearest location on the nearest edge, only the first of them is blended into the network. This is for the same reasons as explained before: in the network structure there is no clear approach for dealing with duplicated nodes. By arranging your table of POIs with dplyr::arrange() before blending you can influence which (type of) POI is given priority in such cases. There is also the option to set ignore_duplicates = FALSE. Then, duplicated projections will result in duplicated nodes, but these duplicates are isolated (i.e. not connected to the rest of the network). You can then use the morpher to_spatial_unique() to merge spatially duplicated nodes into one, specifying how their attributes should be combined. See here for an example.
Secondly: when a single POI has multiple nearest edges, it is only blended into the first of these edges. Therefore, it might be a good idea to run the to_spatial_subdivision() morpher after blending, such that intersecting but unconnected edges get connected.
Lastly: it is important to be aware of floating point precision. See the discussion in this GitHub issue for more background. In short: due to internal rounding of rational numbers in R it is actually possible that even the intersection point between two lines is not evaluated as intersecting those lines themselves. Sounds confusing? It is! But see the example below:
# Create two intersecting lines.
p1 = st_point(c(0.53236, 1.95377))
p2 = st_point(c(0.53209, 1.95328))
l1 = st_sfc(st_linestring(c(p1, p2)))
p3 = st_point(c(0.53209, 1.95345))
p4 = st_point(c(0.53245, 1.95345))
l2 = st_sfc(st_linestring(c(p3, p4)))
# The two lines share an intersection point.
st_intersection(l1, l2)
#> Geometry set for 1 feature
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 0.5321837 ymin: 1.95345 xmax: 0.5321837 ymax: 1.95345
#> CRS: NA
#> POINT (0.5321837 1.95345)
# But this intersection point does not intersects the line itself!
st_intersects(l1, st_intersection(l1, l2), sparse = FALSE)
#> [,1]
#> [1,] FALSE
# The intersection point is instead located a tiny bit next to the line.
st_distance(l1, st_intersection(l1, l2))
#> [,1]
#> [1,] 4.310191e-17That is: you would expect an intersection with an edge to be blended into the network even if you set tolerance = 0, but in fact that will not always happen. To avoid having these problems, you can better set the tolerance to a very small number instead of zero.
net = as_sfnetwork(l1)
p = st_intersection(l1, l2)
plot(l1)
plot(l2, col = "grey", lwd = 2, add = TRUE)
plot(st_network_blend(net, p, tolerance = 0), lwd = 2, cex = 2, add = TRUE)
#> Warning: [1m[22m`st_network_blend()` did not blend any points into the network.
#> [36mℹ[39m Increase `tolerance` for a higher tolerance distance.
plot(l1)
plot(l2, col = "grey", lwd = 2, add = TRUE)
plot(st_network_blend(net, p, tolerance = 1e-10), lwd = 2, cex = 2, add = TRUE)
Joining two networks
In the examples above it was all about joining information from external features into a network. But how about joining two networks? This is what the st_network_join() function is for. It takes two sfnetworks as input and makes a spatial full join on the geometries of the nodes data, based on the equals spatial predicate. That means, all nodes from network x and all nodes from network y are present in the joined network, but if there were nodes in x with equal geometries to nodes in y, these nodes become a single node in the joined network. Edge data are combined using a dplyr::bind_rows() semantic, meaning that data are matched by column name and values are filled with NA if missing in either of the networks. The from and to columns in the edge data are updated automatically such that they correctly match the new node indices of the joined network. There is no spatial join performed on the edges. Hence, if there is an edge in x with an equal geometry to an edge in y, they remain separate edges in the joined network.
# Create two networks.
# Create a network.
n1 = st_point(c(0, 0))
n2 = st_point(c(1, 0))
n3 = st_point(c(1, 1))
n4 = st_point(c(0, 1))
e1 = st_sfc(st_linestring(c(n1, n2)))
e2 = st_sfc(st_linestring(c(n2, n3)))
e3 = st_sfc(st_linestring(c(n3, n4)))
neta = st_sf(geometry = c(e1, e2)) |>
as_sfnetwork()
netb = st_sf(geometry = c(e2, e3)) |>
as_sfnetwork()
# Join them into a single network.
joined = st_network_join(neta, netb)
joined#> # A sfnetwork: 4 nodes and 4 edges
#> #
#> # A directed acyclic multigraph with 1 component and spatially explicit edges
#> #
#> # Dimension: XY
#> # Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1
#> # CRS: NA
#> #
#> # Node data: 4 × 1 (active)
#> geometry
#> <POINT>
#> 1 (0 0)
#> 2 (1 0)
#> 3 (1 1)
#> 4 (0 1)
#> #
#> # Edge data: 4 × 3
#> from to geometry
#> <int> <int> <LINESTRING>
#> 1 1 2 (0 0, 1 0)
#> 2 2 3 (1 0, 1 1)
#> 3 2 3 (1 0, 1 1)
#> # ℹ 1 more row
plot(neta, col = "skyblue", pch = 15, cex = 2, lwd = 4)
plot(netb, col = "orange", pch = 18, cex = 2, lty = 2, lwd = 4, add = TRUE)
plot(joined, cex = 2, lwd = 4)
