How to split several lines by one polygon and get the distance of each part using sf?

Question

The question "How to split one line by several polygons" already has an answer here. My question is similar but different: How to split several lines by one polygon?

The context: I have a border and routes followed by vehicles. I want to

1. decompose the routes, cutting them everytime they cross a border, and
2. for each part, know what are the distances performed.

Some routes have overlapping segments: they go somewhere and come back using the same segment.

I would like to stay with sf for this problem, to limit the number of packages I'm using in my project.

** Input data**: Here is a simple example with a square as the border and five routes

library(sf)

# Defining the "border", a polygon
four_corners <- matrix(c(0, 0, 4, 0, 4, 4, 0, 4, 0, 0), 
                       ncol=2, byrow=TRUE)
pts <- list(four_corners)
my_polygon <- st_polygon(pts)
my_polygon_as_a_linestring <- st_cast(my_polygon, "LINESTRING", group_or_split=FALSE)

# Defining the first "route", with 1 intersection
two_points_of_linestring1 <- matrix(c(0, -1, 2, 1),
                                    ncol=2, byrow=TRUE)
my_linestring1 <- st_linestring(two_points_of_linestring1)

# Defining the second "route", with 2 intersections
two_points_of_linestring2 <- matrix(c(1, -1, 5, 3),
                                    ncol=2, byrow=TRUE)
my_linestring2 <- st_linestring(two_points_of_linestring2)

# Defining the third "route", internal but with overlapping segments
five_internal_points <- matrix(c(1, 2, 2, 2, 2, 3, 2, 2, 3, 2),
                      ncol=2, byrow=TRUE)
my_linestring3 <- st_linestring(five_internal_points)

# Defining the fourth "route", with 3 intersections
four_points <- matrix(c(0.5, 0.5, -0.5, 1.5, 0.5, 2.5, -0.5, 3.5),
                      ncol=2, byrow=TRUE)
my_linestring4 <- st_linestring(four_points)

# Defining the fifth "route", with 5 intersections
five_points <- matrix(c(0.25, 3.5, 1.25, 4.5, 2.25, 3.5, 3.25, 4.5, 4.25, 3.5),
                      ncol=2, byrow=TRUE)
my_linestring5 <- st_linestring(five_points)

# Group the five routes in a spatial dataframe with ID
g <- st_sfc(my_linestring1, my_linestring2, my_linestring3, my_linestring4, my_linestring5)
df_routes <- data.frame(id_column = 1:5)
st_geometry(df_routes) <- g

The geodataframe is thus:

  id_column                       geometry
1         1         LINESTRING (0 -1, 2 1)
2         2         LINESTRING (1 -1, 5 3)
3         3 LINESTRING (1 2, 2 2, 2 3, ...
4         4 LINESTRING (0.5 0.5, -0.5 1...
5         5 LINESTRING (0.25 3.5, 1.25 ...

The expected output is:

  id_column distPart1 distPart2 distPart3 distPart4 distPart5 distPart6 distPart7 distPart8
1         1    1.4...    1.4...        NA        NA        NA        NA        NA        NA       
2         2    1.4...    2.8...    1.4...        NA        NA        NA        NA        NA         
3         3         4        NA        NA        NA        NA        NA        NA        NA
4         4    0.7...    0.7...    0.7...    0.7...    0.7...    0.7...        NA        NA 
5         5    0.7...    0.7...    0.7...    0.7...    0.7...    0.7...    0.7...    0.7... 

By using this approach, we decompose all routes by all intersections, which is not correct.

Using st_split does not work with overlapping segments, as described here.

How can I decompose each route by its own intersections?

Answer 1

I found a solution to quantify the length of the overlaping lines as st_crop or st_intersection seems to get rid of those segments (so the final lengths of the segments doesn't match).

In summary, here is the process:

1. Generate a layer of points with each of the nodes of the lines and each of the intersection points with the polygon.
2. Split the lines individualy using those points. The lines will be quantify and merged at the end by their attributes.
3. Add 'IN' our 'OUT' to each of the parts.
4. Merge all the segments again and quantify the length and the proportions (just to notice quickly if the length of the parts match the initial length).

Here is the code (just execute after yours):

# ------------------------------------------------------------------------------
library(lwgeom)
library(mapview)
library(dplyr)

# UPDATE
# SET layers
ls <- st_as_sf(df_routes) %>% mutate(TOTAL=st_length(df_routes))
p <- my_polygon
pl <- my_polygon_as_a_linestring

# ------------------------------------------------------------------------------
# Extract by breaking points
ps1 <- st_cast(ls, 'POINT') # own nodes
ps2 <- st_intersection(ls,pl) # nodes from intersection
ps <- rbind(ps1,ps2)

mapview(list(ls,ps,pl),
        col.regions=list("blue","green",'red'),
        col=list("blue","green",'red')) 

parts <- st_split(ls, ps)
lscut = st_collection_extract(parts,"LINESTRING")

# ------------------------------------------------------------------------------
# IN LINES
# Crop, merge and get total length. Also add zone attribute
lsin <- st_intersection(lscut,p)
lsin <- lsin %>% mutate(LENGTH = st_length(lsin)) %>% mutate(ZONE = 'IN')

# Remove 0's (starting points from outer line)
lsin <- lsin %>% filter(!LENGTH==0)

# Collapse lines to get lengths
lsin <- lsin %>%
        group_by(id_column, ZONE, TOTAL) %>%
        summarize(LENGTH = sum(LENGTH))

# OUT LINES
# Crop, merge and get total length. Also add zone attribute
lsout <- st_difference(lscut,p)
lsout <- lsout %>% mutate(LENGTH = st_length(lsout)) %>% mutate(ZONE = 'OUT')

# Remove 0's (starting points from outer line)
lsout <- lsout %>% filter(!LENGTH==0)

# Collapse lines to get lengths
lsout <- lsout %>%
        group_by(id_column, ZONE, TOTAL) %>%
        summarize(LENGTH = sum(LENGTH))

# Merge all in a new layer
nl <- rbind(lsin,lsout)
nl <- nl %>% arrange(id_column) %>% mutate(PROP=round((LENGTH/TOTAL)*100)) %>% 
        select(id_column, TOTAL, LENGTH, ZONE, PROP)
nl

mapview(list(ls, nl,pl),
        alpha=list(1, 0.2, 0.5),
        col.regions=list("blue","green","red"),
        col=list("blue","green","red"))