CODEWITHSUNDEEP
rspatialclip
Sébastien Wouters
Sébastien Wouters

Reputation: 131

Complex clipping (spatial intersection ?) of polygons and lines in R

I would like to clip (or maybe the right formulation is performing spatial intersection) polygons and lines using a polygon rather than a rectangle, like so: enter image description here

Here is some code to make the polygons for reproducibility and examples:

p1 <- data.frame(x = c(-0.81, -0.45, -0.04, 0.32, 0.47, 0.86, 0.08, -0.46, -1, -0.76), 
                 y = c(0.46, 1, 0.64, 0.99, -0.04, -0.14, -0.84, -0.24, -0.44, 0.12))

p2 <- data.frame(x = c(-0.63, -0.45, -0.2, -0.38, -0.26, -0.82, -0.57, -0.76), 
                 y = c(-0.1, 0.15, -0.17, -0.79, -1, -0.97, -0.7, -0.61))

l1 <- data.frame(x = c(0.1, 0.28, 0.29, 0.52, 0.51, 0.9, 1), 
                 y = c(0.19, -0.15, 0.25, 0.28, 0.64, 0.9, 0.47))

plot.new()
plot.window(xlim = c(-1, 1), ylim = c(-1,1))

polygon(p2$x, p2$y, col = "blue")
polygon(p1$x, p1$y)

lines(l1$x, l1$y)

Upvotes: 0

Views: 315

Answers (1)

Ege Rubak
Ege Rubak

Reputation: 4507

You could use the spatstat package for this. Below the original example is worked through. In spatstat polygons are used as “observation windows” of point patterns, so they are of class owin. It is possible to do set intersection, union etc. with owin objects.

p1 <- data.frame(x = c(-0.81, -0.45, -0.04, 0.32, 0.47, 0.86, 0.08, -0.46, -1, -0.76), 
                 y = c(0.46, 1, 0.64, 0.99, -0.04, -0.14, -0.84, -0.24, -0.44, 0.12))

p2 <- data.frame(x = c(-0.63, -0.45, -0.2, -0.38, -0.26, -0.82, -0.57, -0.76), 
                 y = c(-0.1, 0.15, -0.17, -0.79, -1, -0.97, -0.7, -0.61))

l1 <- data.frame(x = c(0.1, 0.28, 0.29, 0.52, 0.51, 0.9, 1), 
                 y = c(0.19, -0.15, 0.25, 0.28, 0.64, 0.9, 0.47))

In spatstat polygons must be traversed anti-clockwise, so:

library(spatstat)
p1rev <- lapply(p1, rev)
p2rev <- lapply(p2, rev)

W1 <- owin(poly = p1rev)
W2 <- owin(poly = p2rev)
L1 <- psp(x0 = l1$x[-nrow(l1)], y0 = l1$y[-nrow(l1)],
          x1 = l1$x[-1], y1 = l1$y[-1], window = boundingbox(l1))

plot(boundingbox(W1,W2,L1), type= "n", main = "Original")
plot(W2, col = "blue", add = TRUE)
plot(W1, add = TRUE)
plot(L1, add = TRUE)


W2clip <- W2[W1]
L1clip <- L1[W1]

plot(W1, main = "Clipped")
plot(W2clip, col = "blue", add = TRUE)
plot(L1clip, add = TRUE)

Upvotes: 1

Related Questions