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

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