9

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

4 Answers 4

11

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
  • 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.
    – Antonin
    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 Mar 30 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]. Mar 30 at 18:38
4

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

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)

1
1

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)

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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