terra icon indicating copy to clipboard operation
terra copied to clipboard

Split at intersection for lines?

Open MTueting opened this issue 1 year ago • 1 comments

The split function allows us to split a polygon by a line. Currently, the first input needs to be a polygon, and the second needs to be a line. It would be nice to include the option to split a line by another object (line or point).

Something along these lines:

library(sf)

ls1 = st_sfc(st_linestring(matrix(c(0,1,0,1), ncol = 2)))
ls2 = st_sfc(st_linestring(matrix(c(0,1,1,0), ncol = 2)))
ls = c(ls1, ls2)
ls = st_as_sf(ls)
ls$v = c("v1", "v2")
ls$col = c("blue", "red")
ls
plot(st_geometry(ls), axes = T, col = 3:4, lwd = 3)
p <- st_collection_extract(st_intersection(ls), "POINT")
q <- lwgeom::st_split(ls, p)
q <- st_collection_extract(q, "LINESTRING")

plot(st_geometry(q), axes = T, col = 1:4, lwd = 3)

MTueting avatar Dec 15 '23 00:12 MTueting

For someone interested: There is a workaround using terra's ''aggregate'' function and lwgeom's ''st_split''.

library(terra)

line1 <- vect("LINESTRING (0 0, 2 2)")
line2 <- vect("LINESTRING (0 2, 2 0)")
line3 <- vect("LINESTRING (1 0, 1 2)")

# Combine lines into a single SpatVector
lines <- rbind(line1, line2, line3)
lines$id <- 1:nrow(lines)
plot(lines, y = "id", lwd = 2)

# We make this into a function
split_lines <- function(lines) {
  # Aggregate the input lines
  agg_lines <- aggregate(lines)
  # Find the points that were added
  points_before <- as.points(lines)
  points_after <- as.points(agg_lines)
  points_added <- erase(points_after, points_before)
  # Split the lines at the added points
  split_lines <- sf::st_collection_extract(lwgeom::st_split(sf::st_as_sf(lines), sf::st_as_sf(points_added)), "LINESTRING")
  split_lines$new_id <- 1:nrow(split_lines)
  return(vect(split_lines))
}

# We test the function
out = split_lines(lines)
plot(out, y = "new_id", lwd = 2)

MTueting avatar Jul 11 '24 14:07 MTueting

Thanks, now implemented. See:

library(terra)
v <- vect(system.file("ex/lux.shp", package="terra"))
v <- as.lines(v[c(1,3,4)])
line <- vect(matrix(c(5.79, 6.22, 5.75, 6.1, 5.8, 50.14, 50.05, 49.88, 49.85, 49.71), ncol=2), "line")
s <- split(v, line)
plot(s, col=sample(rainbow(35)), lwd=3)
lines(line)

Image

rhijmans avatar Mar 03 '25 02:03 rhijmans