How do I restrict a SpatialPolygonsDataFrame
so that none of the polygons lie outside of an "underlying" SpatialPolygonsDataFrame
?
The application I have in mind involves Voronoi Polygons a la Carson Farmer -- the result of such a decomposition is (typically) by default rectangular, but often it's more visually appealing to conform the result to an underlying geometry (city/county/country border, etc.). Another common application is to delete Alaska and Hawaii from US-wide geometries when only the lower 48 are of interest.
This is similar to this question, but more complicated, since both input geometries consist of multiple polygons -- a simple gIntersection
will result in an unnecessary dissolution of the main geometry of interest.
It's also similar to this question but I couldn't follow exactly what was going on there because there was no reproducible example. I think the following exercise likely subsumes this question.
Here's a simple example to work with:
library(sp)
library(rgeos)
library(maptools)
points<-list("A" = c(0, 0),"B" = c(1/2, 0),"C" = c(1, 0),
"D" = c(0, 1/2), "E" = c(1/2, 1/2), "F" = c(1, 1/2),
"G" = c(0, 1), "H" = c(1/2, 1), "I" = c(1, 1),
"J" = c(0, -1/2), "K" = c(3/2, -1/2), "L" = c(7/24, 3/8),
"M" = c(5/8, 3/8), "N" = c(3/8, 5/8),
"O" = c(5/8, 17/24), "P" = c(3/2, 1))
pts_poly <- function(x)Reduce(rbind, lapply(x, function(y) unlist(points[y])))
poly1_list<-
list("poly11" = list(ID = "A", coords = pts_poly(c("A", "B", "E", "D" ))),
"poly12" = list(ID = "B", coords = pts_poly(c("B", "C", "F", "E"))),
"poly13" = list(ID = "C", coords = pts_poly(c("E", "F", "I", "H"))),
"poly14" = list(ID = "D", coords = pts_poly(c("D", "E", "H", "G"))))
poly2_list<-
list("poly11" = list(ID = "a", coords = pts_poly(c("J", "M", "L"))),
"poly12" = list(ID = "b", coords = pts_poly(c("J", "K", "M"))),
"poly13" = list(ID = "c", coords = pts_poly(c("M", "K", "P"))),
"poly14" = list(ID = "d", coords = pts_poly(c("M", "P", "O"))),
"poly15" = list(ID = "e", coords = pts_poly(c("L", "M", "O", "N"))))
poly1<-SpatialPolygons(lapply(
poly1_list, function(pl){
with(pl, Polygons(list(Polygon(coords)), ID = ID))}))
df1 <- data.frame(var1 = rnorm(length(poly1)), var2 = rnorm(length(poly1)),
row.names = sapply(poly1@polygons, function(x) x@ID))
polydf1 <- SpatialPolygonsDataFrame(poly1, df1)
poly2<-SpatialPolygons(lapply(
poly2_list, function(pl){
with(pl, Polygons(list(Polygon(coords)), ID = ID))}))
df2 <- data.frame(var1 = rnorm(length(poly2)), var2 = rnorm(length(poly2)),
row.names = sapply(poly2@polygons, function(x) x@ID))
polydf2 <- SpatialPolygonsDataFrame(poly2, df2)
Lots of code. Here's what we have visually:
(code for plots:)
par(mfrow = c(1, 3))
par(mar=c(0,0,0,0))
par(oma=c(0,0,0,0))
plot(polydf1, xlim = c(0, 1.1), ylim = c(0, 1.1))
title("Underlying Polygons", line = -15, cex.main=1.8)
text(Reduce(rbind, points[1:9]) + .03,
labels = names(points)[1:9], col = "red")
text(matrix(c(.25, .25, .75, .25,
.75, .75, .25, .75),
ncol = 2, byrow = TRUE), cex = 2,
labels=sapply(polydf1@polygons, slot, name = "ID"))
plot(polydf2, xlim = c(0, 1.6), ylim = c(-.5, 1.1))
title("Polygons to Be Clipped", line = -15, cex.main=1.8)
text(Reduce(rbind, points[-(1:9)]) + .03,
labels = names(points)[-(1:9)], col = "red")
text(coordinates(polydf2), cex = 2,
labels=sapply(polydf2@polygons, slot, name = "ID"))
plot(polydf2)
plot(polydf1, add = TRUE)
title("Overlay of Polygons", line = -15, cex.main=1.8)
What I want is to eliminate the parts of triangles E, F, G and H which lie outside of the the the first set of polygons, A, B, C, and D.
A simple intersection is incorrect:
plot(gIntersection(polydf1, polydf2, byid = TRUE),
main = "Simple Intersection: Overkill")
Here's what my desired output looks like instead (basically eliminating the borders of A, B, C, and D):