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

Question

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)
plot(my_linestring1, add=TRUE)

# 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)
plot(my_linestring2, add=TRUE)

# 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)
plot(my_linestring3, add=TRUE)

# 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?

Answer 1

Here is how you might do that with terra

Your example data

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

Answer 2

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)).