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

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?

• Instead of `my_linestring1 <- st_linestring(three_points)`, it shouldn't be `my_linestring1 <- st_linestring(two_points_of_linestring1)` as `three_points` does not exist? May 20, 2021 at 0:24
• Yes, indeed! Thank you for spotting the mistake and letting me know. I've done the correction in the question. May 20, 2021 at 7:43
• I think that the line 3 is not well planned as it lacks some segments (for example up to down) in order to be a real connected line. I'm close to the solution but writting that line in a correct way It could work... May 20, 2021 at 11:30
• I'm not sure that I understand your comment. The route goes from left to right ((1,2) -> (2,2)), then up ((2, 2) -> (2, 3)), then down ((2, 3) -> (2, 2)), then right again ((2, 2) -> (3, 2)). If I plot it, I see what on the figure in the question. There are no error message. But of course, there is an overlapping segment (which is the goal). May 20, 2021 at 12:32
• It seems that sf packages interprets those overlaping segments as topology errors: gis.stackexchange.com/questions/276419/… We must find a way to avoid it... May 20, 2021 at 14:09

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"))
`````` 