6

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.

Example

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

7

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)

enter image description here

3

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)

plot

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.