2

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:

enter image description here

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?

1
  • 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.
    – Antonin
    Commented May 7, 2021 at 13:57

2 Answers 2

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
2
  • 1
    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.
    – Antonin
    Commented May 5, 2021 at 14:02
  • 2
    See my expanded answer. Commented May 5, 2021 at 15:09
1

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])
}

route with overlapping

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.