# Split a line when it crosses a polygon in R

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

## 4 Answers

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

• How could we extend your approach when there are two or more lines? I've naively tried and the issue is that st_difference is using all the intersection points with all the lines (and not intersection points of line A with line A and intersection points of line B with line B). See gis.stackexchange.com/questions/395751/… for an example. Commented May 6, 2021 at 15:38
• This answers a seperate question i had very well. how do i then subset/select each part of the line, for example only the blue section Commented Mar 30, 2023 at 7:38
• @Salmosalar Use `st_cast()` to convert the `"MULTILINESTRING"` feature to its constituent `"LINESTRING"` parts, and then select the piece(s) you want, like this: `st_cast(C, "LINESTRING")[3]`. Commented Mar 30, 2023 at 18:38

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

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)

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