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:


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

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

  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_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_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:

input polygons

(code for plots:)

par(mfrow = c(1, 3))
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(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")

basic intersection

Here's what my desired output looks like instead (basically eliminating the borders of A, B, C, and D):

desired output


The simplest approach to do this is

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.

| improve this answer | |
  • 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:

enter image description here

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)
| improve this answer | |
  • 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

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.