How to merge groups of touching lines to multilinestring

Question

I would like to join groups of lines into multilinestrings if they touch. In the plot below, I should get two objects of sf class MULTILINESTRING from the 6 lines. What is an easy way of doing this (when having 2 mio lines)? Is there maybe a pre-defined command that let's one do this?

data:

library(sf)
lines <- st_sfc(list(st_linestring(rbind(c(0,0), c(0.5, 0.7))), 
                 st_linestring(rbind(c(0.5, 0.7), c(1, 1))), 
                 st_linestring(rbind(c(2,0), c(1,1))), 
                 st_linestring(rbind(c(1,1), c(1, 3))),
                 st_linestring(rbind(c(3,2), c(4,2))),
                 st_linestring(rbind(c(3,2), c(5,3)))))

I tried the following

library(igraph)
components(graph.adjlist(st_touches(lines)))

graph.adjlist transforms the list of elements that touch each other to a graph object and the function components identifies the components of each cluster.

But is there a way to do connect the multiple LINESTRING into MULTILINESTRING if they touch?

Answer 1

Use tapply over the groups defined by the connectivity to make multilinestrings:

> com = components(graph.adjlist(st_touches(lines)))

> mlines = do.call(st_sfc,
              tapply(lines,
                     com$membership,
                     function(i){
                      st_multilinestring(lines[i])}))

Plotting the two features:

> plot(mlines,col=c("red","blue"),lwd=4)

Answer 2

A Tidyverse alternative:

library("dplyr")  

com = components(graph.adjlist(st_touches(lines)))

# generate a simple features collection and add the group number
lines = st_sf(geometry = lines)
lines$group = com$membership 

# summarise by group
multi_lines = lines |> group_by(group) |>
  summarise()

multi_lines
#> Simple feature collection with 2 features and 1 field
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 0 ymin: 0 xmax: 5 ymax: 3
#> CRS:           NA
#> # A tibble: 2 × 2
#>   group                                                 geometry
#>   <dbl>                                        <MULTILINESTRING>
#> 1     1 ((0 0, 0.5 0.7), (0.5 0.7, 1 1), (2 0, 1 1), (1 1, 1 3))
#> 2     2                                 ((3 2, 4 2), (3 2, 5 3))

plot(multi_lines$geometry, col = multi_lines$group)

Although it seems that aggregate would also be a possibility: aggregate(lines, by = list(lines$group), FUN = unique)

^{Created on 2023-08-11 with reprex v2.0.2}