# Clipping a SpatialPoygonsDataFrame to another SpatialPolygonsDataFrame

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)
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): The simplest approach to do this is

``````library(raster)
x <- crop(polydf2, polydf1)
``````

This will do the geometric intersection, but also assure the attributes are in order. See `?'raster-package'` (section XIV) for more functions that deal with polygon overlay.

• thanks! accepting my own answer only because I'm a package minimalist – MichaelChirico Nov 15 '15 at 8:28

The answer is pretty simple once you see it, but I spent half a day banging my head on the desk before realizing the simple adjustment that we need to make; I couldn't find it anywhere else online, so I wanted to leave this here for the reference of anyone else struggling like myself.

The key is to eliminate the borders of the first set of polygons prior to intersecting with the important polygons:

``````clipped_polys <-
gIntersection(polydf2, gUnaryUnion(polydf1),
byid = TRUE, id = sapply(polydf2@polygons, slot, name = "ID"))
``````

To make sure the polygons retain their original names, we also specify the ID slot ourselves. Here's the output: And of course we can add back any data associated with `poly2` readily as pointed out in this answer:

``````clipped_polys <- SpatialPolygonsDataFrame(clipped_polys, polydf2@data)
``````
• This works for this case, but it is not a good generic solution. For example, try `polydf1 <- polydf1[3:4, ]` etc. The last line will fail. – Robert Hijmans Nov 15 '15 at 18:39