Split a line when it crosses a polygon in R

Question

This should be a simple problem but I can't seem to find a solution.

I have a line which crosses several polygons and I want to split it into multiple lines at the intersection of each polygon using R. In ArcGIS the Union tool can achieve this but Union in R seems to work differently. I using sf but could also use sp if needed.

The image shows a single line (red), which crosses the two polygons and a gap between the polygons. The result should be three lines (red,green,blue).

Answer 1

Here is an answer that applies sf package functions to the reproducible data kindly provided by rcs.

library(sf)
A <- st_as_sfc("LINESTRING(458.1 768.23, 455.3 413.29, 522.3 325.77, 664.8 282.01, 726.3 121.56)")
B <- st_as_sfc("MULTIPOLYGON(((402.2 893.03, 664.8 800.65, 611.7 666.13, 368.7 623.99, 215.1 692.06, 402.2 893.03)), ((703.9 366.29, 821.2 244.73, 796.1  25.93, 500.0 137.76, 703.9 366.29)))")

## Convert the MULTIPOLYGON to a MULTILINESTRING object
BB <- st_cast(B, "MULTILINESTRING", group_or_split=FALSE)

## Break LINESTRING A into segments by using:
## - st_intersection() to find points at which lines features intersect
## - st_buffer() to convert points to tiny polygons with some width
## - st_difference() to break line up into sections not overlapping tiny polygons
C <- st_difference(A, st_buffer(st_intersection(A,BB), dist=1e-12))

## Check that it works:
plot(B, col="grey")
plot(st_cast(C, "LINESTRING"),
     col=c("red2", "springgreen3", "dodgerblue"), lwd=3, add=TRUE)

Answer 2

Here is a sp/rgeos solution, mainly based on the linear referencing functions gInterpolate and gProject:

Reproducible Example:

library("sp")
library("rgeos")

# Return a linestring being a substring of the first argument, starting and
# ending at the given fractions of total 2d length. Second and third arguments
# are numeric values between 0 and 1.
line_substring <- function(spgeom, pos1=0, pos2=1) {

  stopifnot(inherits(spgeom, "SpatialLines") ||
            inherits(spgeom, "SpatialLinesDataFrame"))

  val_line <- gProject(spgeom, as(spgeom, "SpatialPoints"), normalized=TRUE)
  ind <- (val_line >= pos1) & (val_line <= pos2)

  res <- list(gInterpolate(spgeom, pos1, normalized=TRUE),
              as(spgeom, "SpatialPoints")[ind, ],
              gInterpolate(spgeom, pos2, normalized=TRUE))
  as(do.call(rbind, res), "SpatialLines")
}


# example data
l <- readWKT("LINESTRING(458.1 768.23, 455.3 413.29, 522.3 325.77, 664.8 282.01, 726.3 121.56)")
p <- readWKT("MULTIPOLYGON(((402.2 893.03, 664.8 800.65, 611.7 666.13, 368.7 623.99, 215.1 692.06, 402.2 893.03)), ((703.9 366.29, 821.2 244.73, 796.1  25.93, 500.0 137.76, 703.9 366.29)))")

# get intersection points
p_intersect <- gIntersection(as(p, "SpatialLines"), l)

# project intersection points to line
line_dist <- gProject(l, p_intersect, normalized=TRUE)
line_dist <- c(0, line_dist, 1)

# list `res` contains the resulting lines
res <- list()
for (i in seq_len(length(line_dist)-1)) {
  res[[i]] <- line_substring(l, line_dist[i], line_dist[i+1])
}


plot(p, axes=TRUE, border="gray")
for (i in seq_along(res)) plot(res[[i]], col=i+1, lwd=4, add=TRUE)
plot(l, add=TRUE, lty=2)

Answer 3

Already an older question, but for those who read this and are interested: you can also use the package lwgeom for this. Written by the author of sf, it gives access to selected functions found in liblwgeom, the light-weight geometry library used by PostGIS.

It contains a function st_split, which will split the line at the intersections and returns it as a geometrycollection. Extracting the elements from this geometrycollection gives you the separate parts of the splitted line:

library(sf)
#> Linking to GEOS 3.7.1, GDAL 2.2.2, PROJ 4.9.2
library(lwgeom)
#> Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.7.1, PROJ 4.9.2

A <- st_as_sfc("LINESTRING(458.1 768.23, 455.3 413.29, 522.3 325.77, 664.8 282.01, 726.3 121.56)")
B <- st_as_sfc("MULTIPOLYGON(((402.2 893.03, 664.8 800.65, 611.7 666.13, 368.7 623.99, 215.1 692.06, 402.2 893.03)), ((703.9 366.29, 821.2 244.73, 796.1  25.93, 500.0 137.76, 703.9 366.29)))")

A_new = st_split(A, B)
A_new
#> Geometry set for 1 feature 
#> geometry type:  GEOMETRYCOLLECTION
#> dimension:      XY
#> bbox:           xmin: 455.3 ymin: 121.56 xmax: 726.3 ymax: 768.23
#> CRS:            NA
#> GEOMETRYCOLLECTION (LINESTRING (458.1 768.23, 4...

A_new = st_collection_extract(A_new, "LINESTRING")
A_new
#> Geometry set for 3 features 
#> geometry type:  LINESTRING
#> dimension:      XY
#> bbox:           xmin: 455.3 ymin: 121.56 xmax: 726.3 ymax: 768.23
#> CRS:            NA
#> LINESTRING (458.1 768.23, 457.083 639.317)
#> LINESTRING (457.083 639.317, 455.3 413.29, 522....
#> LINESTRING (636.4665 290.7109, 664.8 282.01, 72...

colors = c("red", "green", "blue")
plot(B)
for (i in seq_along(A_new)) {
  plot(A_new[i], col = colors[i], add = T)
}

^{Created on 2020-10-28 by the reprex package (v0.3.0)}

Answer 4

You can do this with terra like this:

library(terra)
A <- vect("LINESTRING(458.1 768.23, 455.3 413.29, 522.3 325.77, 664.8 282.01, 726.3 121.56)")
B <- vect("MULTIPOLYGON(((402.2 893.03, 664.8 800.65, 611.7 666.13, 368.7 623.99, 215.1 692.06, 402.2 893.03)), ((703.9 366.29, 821.2 244.73, 796.1  25.93, 500.0 137.76, 703.9 366.29)))")

x <- disaggregate(intersect(A, B))
y <- erase(A, B)
s <- rbind(x, y)

plot(B)
lines(s, col=c("red", "green", "blue"), lwd=3)