# Split lines by a polygon - Why are "internal" lines split in parts using sf/st_split

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.

Simple reproducible example, part 1: input data: as input, one square as polygon for the border, three routes as linestrings:

``````library(sf)
library(lwgeom)

# Defining the "border", a polygon
four_corners <- matrix(c(0, 0, 1, 0, 1, 1, 0, 1, 0, 0),
ncol=2, byrow=TRUE)
pts <- list(four_corners)
my_polygon <- st_polygon(pts)
plot(my_polygon)

# Defining the first "route", a linestring
three_points <- matrix(c(-1, 0, 0.5, 0.5, 2, 0),
ncol=2, byrow=TRUE)
my_linestring1 <- st_linestring(three_points)
print(my_linestring1)

# Defining the second "route", a linestring
two_points <- matrix(c(0.5, -0.5, 0.5, 0.25),
ncol=2, byrow=TRUE)
my_linestring2 <- st_linestring(two_points)
print(my_linestring2)

# Defining the third "route", a linestring
four_points <- matrix(c(0.25, 0.75, 0.25, 0.6, 0.75, 0.6, 0.75, 0.75),
ncol=2, byrow=TRUE)
my_linestring3 <- st_linestring(four_points)
print(my_linestring3)

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

`df_routes` looks like this:

``````  id_column                       geometry
1         1 LINESTRING (-1 0, 0.5 0.5, ...
2         2 LINESTRING (0.5 -0.5, 0.5 0...
3         3 LINESTRING (0.25 0.75, 0.25...
``````

With one internal route, one route crossing the border once, and one route crossing the border twice: Simple reproducible example, part 2: the solution

``````# Split the route in different subroutes cut by the border
splitted_routes <- st_split(df_routes, my_polygon)

# Decompose the GEOMETRYCOLLECTION of LINESTRINGs in a list of LINESTRINGs
splitted_routes_as_single_linestring <- st_collection_extract(splitted_routes, "LINESTRING")

# Get the length of the LINESTRING
list_of_length = st_length(splitted_routes_as_single_linestring)

# Remove the geometry, in order to get a dataframe as output
splitted_routes_as_single_linestring\$geometry <- NULL

# Just merge the length in a column distance in the original dataframe
ID_and_length <- cbind(splitted_routes_as_single_linestring, distance = list_of_length)
``````

This approach works on this simple example and generates the distances of each part on each side of the border:

``````    id_column distance
1           1 1.054093
1.1         1 1.054093
1.2         1 1.054093
2           2 0.500000
2.1         2 0.250000
3           3 0.800000
``````

In particular, the 3rd route is internal and the output is simply one distance corresponding to the complete distance of the original route.

Now, when using this approach on a real life example, some internal routes are decomposed in small pieces. They should be decomposed as one, but somehow, my approach is not working. When I am extracting LINESTRINGS out of the GEOMETRYCOLLECTIONS coming from st_split, I got several objects, which was not my goal/expectation.

What's happening?

• I've decomposed my question in two parts. Here, it is just about st_split splitting routes with overlapping segments inside a polygon. A general question about how to compute the length of segments of a line cut by a polygon is now here: gis.stackexchange.com/questions/396066/…, in particular when there are several routes. May 7, 2021 at 13:57

Here is how you might do that with `terra`

``````library(terra)
c4 <- matrix(c(0, 0, 1, 0, 1, 1, 0, 1, 0, 0), ncol=2, byrow=TRUE)
p <- vect(c4, "polygon")
ln1 <- vect(matrix(c(-1, 0, 0.5, 0.5, 2, 0), ncol=2, byrow=TRUE), "lines")
ln2 <- vect(matrix(c(0.5, -0.5, 0.5, 0.25), ncol=2, byrow=TRUE), "lines")
ln3 <- vect(matrix(c(0.25, 0.75, 0.25, 0.6, 0.75, 0.6, 0.75, 0.75), ncol=2, byrow=TRUE), "lines")
lns <- rbind(ln1, ln2, ln3)
plot(lns)
lines(p, col="blue")
``````

Solution

``````lns\$id <- 1:nrow(lns)
x <- intersect(lns, p)
e <- erase(lns, p)

x\$inside <- TRUE
e\$inside <- FALSE
xe <- rbind(x, e)
xe\$distance <- perimeter(xe)
# in future versions
# xe\$distance <- perim(xe)

values(xe)
#  id inside distance
#1  1   TRUE 1.054093
#2  2   TRUE 0.250000
#3  3   TRUE 0.800000
#4  1  FALSE 2.108185
#5  2  FALSE 0.500000
``````

Based on the comment, to get what you are after

``````lns\$id <- 1:nrow(lns)
x <- intersect(lns, p)
e <- erase(lns, p)
xe <- rbind(x, e)
d <- disaggregate(xe)
d\$distance <- perimeter(d)
values(d)
#  id distance
#1  1 1.054093
#2  2 0.250000
#3  3 0.800000
#4  1 1.054093
#5  1 1.054093
#6  2 0.500000
``````
• Thank you for your answer. I was a bit imprecise the "context" part of my question, writing about "inside" and "outside". In fact, I don't only need distance inside and distance outside, I would like to cut the route everytime it crosses the border. Concretely: my output is of size 6, yours is of size 5. I am sure that your solution with terra can be adapted to do this, but I don't know (yet) how. If you have hint, it would be with pleasure. I'll look on my side. May 5, 2021 at 14:02
• See my expanded answer. May 5, 2021 at 15:09

After a bit of research, I found where the problem is coming from in my real life example: some segments of the linestring are overlapping. When this appears, the linestring is decomposed in several pieces, even when it is indeed in the polygon.

Simple example with overlapping segments

``````library(sf)
library(lwgeom)

# Define the border, a polygon
four_corners <- matrix(c(0, 0, 4, 0, 4, 3, 0, 3, 0, 0),
ncol=2, byrow=TRUE)
pts <- list(four_corners)
my_polygon <- st_polygon(pts)

# Define a route in the middle of the polygon, not crossing it
x_coord = c(1, 2, 2, 2, 3)
y_coord = c(1, 1, 2, 1, 1)
coordinate_matrix <- cbind(x_coord, y_coord)
one_route = st_linestring(coordinate_matrix)

# Split the routes in different subroutes cut by the border (st_split)
splitted_routes <- st_split(one_route, my_polygon)
print(length(splitted_routes))
for (parts in 1:length(splitted_routes)) {
print(splitted_routes[parts])
}
`````` I was expecting length(splitted_routes) to be 1 (`LINESTRING (1 1, 2 1, 2 2, 2 1, 3 1)`) and not 3 (`LINESTRING (1 1, 2 1)`, `LINESTRING (2 1, 2 2)` and `LINESTRING (2 1, 3 1)`).