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.
p1 = st_point(c(4151358, 3208045))
p2 = st_point(c(4151340, 3207120))
p3 = st_point(c(4151856, 3207106))
p4 = st_point(c(4151874, 3208031))
poly = st_multipoint(c(p1, p2, p3, p4)) %>%
st_cast("POLYGON") %>%
st_sfc(crs = 3035)
net = as_sfnetwork(roxel) %>%
st_transform(3035)
filtered = st_filter(net, poly, .pred = st_intersects)
plot(net, col = "grey")
plot(poly, border = "red", lty = 4, lwd = 4, add = TRUE)
plot(net, col = "grey")
plot(filtered, add = TRUE)
filtered = net %>%
activate("edges") %>%
st_filter(poly, .pred = st_intersects)
plot(net, col = "grey")
plot(poly, border = "red", lty = 4, lwd = 4, add = TRUE)
plot(net, col = "grey")
plot(filtered, add = TRUE)
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.
filtered = net %>%
activate("edges") %>%
st_filter(poly, .pred = st_intersects) %>%
activate("nodes") %>%
filter(!node_is_isolated())
plot(net, col = "grey")
plot(poly, border = "red", lty = 4, lwd = 4, add = TRUE)
plot(net, col = "grey")
plot(filtered, add = TRUE)
Filtering can also be done with other predicates.
point = st_centroid(st_combine(net))
filtered = net %>%
activate("nodes") %>%
st_filter(point, .predicate = st_is_within_distance, dist = 500)
plot(net, col = "grey")
plot(point, col = "red", cex = 3, pch = 20, add = TRUE)
plot(net, col = "grey")
plot(filtered, add = TRUE)
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.
filtered = net %>%
activate("edges") %>%
filter(edge_intersects(poly)) %>%
activate("nodes") %>%
filter(!node_is_isolated())
plot(net, col = "grey")
plot(poly, border = "red", lty = 4, lwd = 4, add = TRUE)
plot(net, col = "grey")
plot(filtered, add = TRUE)
A nice application of this in road networks is to find underpassing and overpassing roads (i.e. edges that cross other edges but are not connected at that point). As we can see in the example below, such roads are not present in our Roxel data, which results in a network without edges.
The tidygraph::.E() function used in the example makes
it possible to directly access the complete edges table inside verbs. In
this case, that means that for each edge we evaluate if it crosses with
any other edge in the network. Similarly, we can use
tidygraph::.N() to access the nodes table and
tidygraph::.G() to access the network object as a
whole.
#> # A sfnetwork with 701 nodes and 0 edges
#> #
#> # CRS: EPSG:3035
#> #
#> # A rooted forest with 701 trees with spatially explicit edges
#> #
#> # Edge data: 0 × 5 (active)
#> # ℹ 5 variables: from <int>, to <int>, name <chr>, type <fct>,
#> # geometry <GEOMETRY [m]>
#> #
#> # Node data: 701 × 1
#> geometry
#> <POINT [m]>
#> 1 (4151491 3207923)
#> 2 (4151474 3207946)
#> 3 (4151398 3207777)
#> # ℹ 698 more rows
If you just want to store the information about the investigated
spatial relation, without filtering the network, you can also use the
spatial node and edge query functions inside a
dplyr::mutate() verb.
#> # A sfnetwork with 701 nodes and 851 edges
#> #
#> # CRS: EPSG:3035
#> #
#> # A directed multigraph with 14 components with spatially explicit edges
#> #
#> # Node data: 701 × 2 (active)
#> geometry in_poly
#> <POINT [m]> <lgl>
#> 1 (4151491 3207923) TRUE
#> 2 (4151474 3207946) TRUE
#> 3 (4151398 3207777) TRUE
#> 4 (4151370 3207673) TRUE
#> 5 (4151408 3207539) TRUE
#> 6 (4151421 3207592) TRUE
#> # ℹ 695 more rows
#> #
#> # Edge data: 851 × 5
#> from to name type geometry
#> <int> <int> <chr> <fct> <LINESTRING [m]>
#> 1 1 2 Havixbecker Strasse residential (4151491 3207923, 4151474 32079…
#> 2 3 4 Pienersallee secondary (4151398 3207777, 4151390 32077…
#> 3 5 6 Schulte-Bernd-Strasse residential (4151408 3207539, 4151417 32075…
#> # ℹ 848 more rows
Besides predicate query functions, you can also use the coordinate query functions for spatial filters on the nodes. For example:
v = 4152000
l = st_linestring(rbind(c(v, st_bbox(net)["ymin"]), c(v, st_bbox(net)["ymax"])))
filtered_by_coords = net %>%
activate("nodes") %>%
filter(node_X() > v)
plot(net, col = "grey")
plot(l, col = "red", lty = 4, lwd = 4, add = TRUE)
plot(net, col = "grey")
plot(filtered_by_coords, add = TRUE)
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(poly) %>%
activate("nodes") %>%
filter(!node_is_isolated())
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
plot(net, col = "grey")
plot(poly, border = "red", lty = 4, lwd = 4, add = TRUE)
plot(net, col = "grey")
plot(clipped, add = TRUE)
Note: Neither of the clipping function currently works well with
undirected networks!
