Split a line when it crosses a polygon in R

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)

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)