Best way to construct clusters of consequtive edges?

[mkvasnicka]

Hello!

I need to start from an edge and find its neighbors, their neighbors, and so on, to some distance. These edges then constitute a cluster (other criteria would be used in practice). Could you hint to me the best way to do it, or more precisely, the best/fastest way to find an edge's neighboring edge? I created a simple function using edge_touches(), but my function is slow.

My approach so far:

My function:

make_cluster <- function(net, core_ids, steps = 1) {
    update_cluster <- function(net, clust_id, cluster_id) {
        net_cluster <- net |> pull(cluster)
        mutate(net, cluster = if_else(id %in% clust_id, cluster_id, net_cluster))
    }

    cluster_id <- min(core_ids)
    net_active <- active(net)
    net <- net |> activate(edges) |>
        update_cluster(core_ids, cluster_id)
    clust <- net |> filter(id == core_ids)
    for (k in seq_len(steps)) {
        clust <- net |>
            filter(edge_touches(clust), !(id %in% core_ids))
        clust_id <- clust |> pull(id)
        net <- update_cluster(net, clust_id, cluster_id)
        core_ids <- unique(c(core_ids, clust_id))
    }
    net #|> activate(net_active)  # not certain how to reactivate
}

A testable code:

library(sf)
library(spNetwork)
library(sfnetworks)
library(tidyverse)
library(tmap)


tmap_mode("view")

roxel <- st_transform(roxel, 31466)
roxel <- lixelize_lines(roxel, 5, 3)
roxel <- mutate(roxel, id = row_number(), cluster = id)
roxy <- as_sfnetwork(roxel)


core_id <- c(3164L)  # 269L)


make_cluster <- function(net, core_ids, steps = 1) {
    update_cluster <- function(net, clust_id, cluster_id) {
        net_cluster <- net |> pull(cluster)
        mutate(net, cluster = if_else(id %in% clust_id, cluster_id, net_cluster))
    }

    cluster_id <- min(core_ids)
    net_active <- active(net)
    net <- net |> activate(edges) |>
        update_cluster(core_ids, cluster_id)
    clust <- net |> filter(id == core_ids)
    for (k in seq_len(steps)) {
        clust <- net |>
            filter(edge_touches(clust), !(id %in% core_ids))
        clust_id <- clust |> pull(id)
        net <- update_cluster(net, clust_id, cluster_id)
        core_ids <- unique(c(core_ids, clust_id))
    }
    net #|> activate(net_active)
}


system.time(
    mm <- make_cluster(roxy, core_id, 60)
)


tm_shape(roxel) + tm_lines(lwd = 3) +
    tm_shape(mm |> activate(edges) |> filter(cluster == core_id) |> st_as_sf()) +
    tm_lines(lwd = 3, col = "green") +
    tm_shape(roxel |> filter(id == core_id)) + tm_lines(col = "red")

Many thanks for any help.
Michal

[agila5]

Hi @mkvasnicka. I think that if you want to subset a cluster of consecutive edges, then you can readapt the following code:

Load packages

library(sf)
library(spNetwork)
library(tidygraph)
library(dplyr)
library(tmap)
library(sfnetworks) # installed from develop branch
library(igraph)

Prepare data

roxel <- st_transform(roxel, 31466)
roxel <- lixelize_lines(roxel, 5, 3)
roxel <- mutate(roxel, id = row_number(), cluster = id)
roxy <- as_sfnetwork(roxel, directed = FALSE)
core_id <- 3164L

Define make_cluster

make_cluster <- function(net, core_ids, steps = 1) {
  update_cluster <- function(net, clust_id, cluster_id) {
    net_cluster <- net |> pull(cluster)
    mutate(net, cluster = if_else(id %in% clust_id, cluster_id, net_cluster))
  }
  
  cluster_id <- min(core_ids)
  net_active <- active(net)
  net <- net |> activate(edges) |>
    update_cluster(core_ids, cluster_id)
  clust <- net |> filter(id == core_ids)
  for (k in seq_len(steps)) {
    clust <- net |>
      filter(edge_touches(clust), !(id %in% core_ids))
    clust_id <- clust |> pull(id)
    net <- update_cluster(net, clust_id, cluster_id)
    core_ids <- unique(c(core_ids, clust_id))
  }
  net #|> activate(net_active)
}

An alternative approach based on ego and line_graph. The idea here is to define a linegraph (see ?make_line_graph) and then use ego (see ?ego) to compute the index of the neighbouring edges.

subset_ego <- function(data, nodes, order) {
  # Transform into linegraph and compute neighbor for "node" and "order"
  idx <- ego(make_line_graph(data), order = order, nodes = nodes)
  
  # Subset edges using idx
  data %E>% slice(as_ids(idx[[1]])) 
}

Compare timing and results

system.time({
  mm <- make_cluster(roxy, core_ids = core_id, steps = 60)
})
#>    user  system elapsed 
#>   25.08    0.17   25.72
system.time({
  mm_ego <- subset_ego(roxy, nodes = core_id, order = 60)
})
#>    user  system elapsed 
#>    0.08    0.00    0.08

mm %E>% filter(cluster == core_id) %E>% st_as_sf()
#> Simple feature collection with 731 features and 7 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 2605229 ymin: 5758468 xmax: 2605762 ymax: 5758961
#> Projected CRS: DHDN / 3-degree Gauss-Kruger zone 2
#> # A tibble: 731 x 8
#>     from    to lineID name         type                   geometry    id cluster
#>    <int> <int>  <int> <chr>        <fct>          <LINESTRING [m]> <int>   <int>
#>  1     9    10      8 Pienersallee seco~ (2605390 5758916, 260539~     8    3164
#>  2    10    11      9 Pienersallee seco~ (2605390 5758911, 260538~     9    3164
#>  3    11    12     10 Pienersallee seco~ (2605389 5758906, 260538~    10    3164
#>  4    12    13     11 Pienersallee seco~ (2605389 5758901, 260538~    11    3164
#>  5    13    14     12 Pienersallee seco~ (2605388 5758896, 260538~    12    3164
#>  6    14    15     13 Pienersallee seco~ (2605388 5758891, 260538~    13    3164
#>  7    15    16     14 Pienersallee seco~ (2605387 5758886, 260538~    14    3164
#>  8    16    17     15 Pienersallee seco~ (2605387 5758881, 260538~    15    3164
#>  9    17    18     16 Pienersallee seco~ (2605386 5758876, 260538~    16    3164
#> 10    18    19     17 Pienersallee seco~ (2605386 5758871, 260538~    17    3164
#> # ... with 721 more rows
mm_ego %E>% st_as_sf()
#> Simple feature collection with 731 features and 7 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 2605229 ymin: 5758468 xmax: 2605762 ymax: 5758961
#> Projected CRS: DHDN / 3-degree Gauss-Kruger zone 2
#> # A tibble: 731 x 8
#>     from    to lineID name         type                   geometry    id cluster
#>    <int> <int>  <int> <chr>        <fct>          <LINESTRING [m]> <int>   <int>
#>  1     9    10      8 Pienersallee seco~ (2605390 5758916, 260539~     8       8
#>  2    10    11      9 Pienersallee seco~ (2605390 5758911, 260538~     9       9
#>  3    11    12     10 Pienersallee seco~ (2605389 5758906, 260538~    10      10
#>  4    12    13     11 Pienersallee seco~ (2605389 5758901, 260538~    11      11
#>  5    13    14     12 Pienersallee seco~ (2605388 5758896, 260538~    12      12
#>  6    14    15     13 Pienersallee seco~ (2605388 5758891, 260538~    13      13
#>  7    15    16     14 Pienersallee seco~ (2605387 5758886, 260538~    14      14
#>  8    16    17     15 Pienersallee seco~ (2605387 5758881, 260538~    15      15
#>  9    17    18     16 Pienersallee seco~ (2605386 5758876, 260538~    16      16
#> 10    18    19     17 Pienersallee seco~ (2605386 5758871, 260538~    17      17
#> # ... with 721 more rows

^{Created on 2022-06-19 by the reprex package (v2.0.1)}

A few notes:

1. I'm not sure what is the expected output, so you may need to tweak that code.
2. As you can see, the code does not update the columns id and cluster (they are ignored for the second approach).
3. If you compute neighbours using the edge_touches() function, then you are implicitly assuming that the input network is undirected. This is an explicit assumption for the second approach (meaning that you will get different results), so I had to set directed = FALSE in the as_sfnetwork() constructor.