How to perform a vector overlay of two SpatialPolygonsDataFrame objects?

Question

I have two GIS layers -- call them Soils and Parcels -- stored as SpatialPolygonsDataFrames (SPDFs), and I would like to "overlay" them, in the sense described here.

The result of the overlay operation should be a new SPDF in which:

1. The SpatialPolygons component contains polygons formed by the intersection of the two layers. (Think all of the atomic polygons formed by overlaying two mylars on an overhead projector).

2. The data.frame component records the attributes of the Soils and Parcels polygons within which each atomic polygon falls.

My question(s): Is there an existing R function that does this? (I would even by happy to learn of a function that just gets the SpatialPolygons component right, calculating the atomic polygons formed from the intersection of the two layers.) I feel like rgeos should have a function that does (1) at least, but it seems not to...

Here is a figure that may help make it clearer what I'm after, followed by code that creates the Soils and Parcels layers shown in the figure.

library(rgeos)

## Just a utility function to construct the example layers.
flattenSpatialPolygons <- function(SP) {
    nm <- deparse(substitute(SP))
    AA <- unlist(lapply(SP@polygons, function(X) X@Polygons))
    SpatialPolygons(lapply(seq_along(AA),
                           function(X) Polygons(AA[X], ID=paste0(nm, X))))
}

## Example Soils layer
Soils <-
local({
    A <- readWKT("MULTIPOLYGON(((3 .5,7 1,7 2,3 1.5,3 0.5), (3 1.5,7 2,7 3,3 2.5,3 1.5)))")
    AA <- flattenSpatialPolygons(A)
    SpatialPolygonsDataFrame(AA,
           data.frame(soilType = paste0("Soil_", LETTERS[seq_along(AA)]),
                      row.names = getSpPPolygonsIDSlots(AA)))
})

## Example Parcels layer
Parcels <-
local({
    B <- readWKT("MULTIPOLYGON(((0 0,2 0,2 3,0 3,0 0),(2 0,4 0,4 3,2 3,2 0)),((4 0,6 0,6 3,4 3,4 0)))")
    BB <- flattenSpatialPolygons(B)
    SpatialPolygonsDataFrame(BB,
           data.frame(soilType = paste0("Parcel_", seq_along(BB)),
                      row.names = getSpPPolygonsIDSlots(BB)))
})

Answer 1

Since January of 2014, the raster package includes a union() function that makes this easy-peasy:

library(raster)
Intersects <- Soils + Parcels  ## Shorthand for union(Soils, Parcels)

## Check that it work
data.frame(Intersects)
## soilType.1 soilType.2
## 1     Soil_A       <NA>
## 2     Soil_B       <NA>
## 3       <NA>   Parcel_1
## 4       <NA>   Parcel_2
## 5       <NA>   Parcel_3
## 6     Soil_A   Parcel_2
## 7     Soil_A   Parcel_3
## 8     Soil_B   Parcel_2
## 9     Soil_B   Parcel_3
plot(Intersects, col = blues9)

Answer 2

Since January 2014, the solution posted here has been completely superseded by the raster::union() function, demonstrated in the officially accepted answer above.

---

This gets the "atomic" polygons formed by overlaying two different SpatialPolygons layers, solving Part 1 of my question.

library(rgeos)

gFragment <- function(X, Y) {
    aa <- gIntersection(X, Y, byid = TRUE)
    bb <- gDifference(X, gUnionCascaded(Y), byid = T)
    cc <- gDifference(Y, gUnionCascaded(X), byid = T)
    ## Note: testing for NULL is needed in case any of aa, bb, or cc is empty,
    ## as when X & Y don't intersect, or when one is fully contained in the other
    SpatialPolygons(c(if(is.null(aa)) NULL else aa@polygons,
                      if(is.null(bb)) NULL else bb@polygons,
                      if(is.null(cc)) NULL else cc@polygons)) 
}

## Try it out
Fragments <- gFragment(Parcels, Soils)
plot(Fragments, col=blues9)

And this extracts the ids (if any) of the polygons in the two input layers overlain by each "atomic" polygon outputted by gFragment() above, solving part 2 of my question.

getAttributes <- function(Fragments, Layer1, Layer2, eps = 0) {
    ## Function to extract attributes from polygon layers
    ## overlain by fragments
    OVER <- function(AA, BB) {
        X <- gRelate(AA, BB, byid = TRUE, pattern="2********")
        ii <- sapply(seq_len(ncol(X)),
                    function(i) {
                        A <- which(X[,i])
                        if(!length(A)) NA else A
                    })
        rownames(X)[ii]
    }
    ## First need to (very slightly) trim Fragments b/c otherwise they
    ## tend to (very slightly) overlap adjacent ring(s)
    Frags <- gBuffer(Fragments, width = -eps, byid = TRUE)
    ## Construct data.frame of attributes
    df <- data.frame(OVER(Frags, Layer1), 
                     OVER(Frags, Layer2),
                     row.names = names(Fragments))
    names(df) <- c(deparse(substitute(Layer1)), deparse(substitute(Layer2)))
    ## Add data.frame to SpatialPolygons object
    SpatialPolygonsDataFrame(Fragments, data=df)
}

FragmentsDF <- getAttributes(Fragments = Fragments,
                             Layer1 = Parcels,
                             Layer2 = Soils)

## A few ways to examine the results
data.frame(FragmentsDF, row.names=NULL)
#   Parcels Soils
# 1      B2    A1
# 2      B2    A2
# 3      B3    A1
# 4      B3    A2
# 5      B1  <NA>
# 6      B2  <NA>
# 7      B3  <NA>
# 8    <NA>    A1
# 9    <NA>    A2
spplot(FragmentsDF, zcol="Soils", col.regions=blues9[3:4])
spplot(FragmentsDF, zcol="Parcels", col.regions=grey.colors(3))

Edit:

Note that this code may fail if any of your input polygons have id's named "1". In that case, one workaround is to rename the id's, perhaps by doing something like Parcels <- spChFIDs(Parcels, paste0("pp", row.names(Parcels))).

Answer 3

Here's my crack, this just gives a list of the pieces with the Parcel (Parcels->Soils) data. You still need to get the attributes from the Soils object, and then do a similar job with the "differences", and then vice versa (Soils->Parcels) so you can have any kinds of overlap relations.

intersects <- list()

## find all intersections (NULLs do nothing to the result)
for (i in 1:nrow(Soils)) {
  for (j in 1:nrow(Parcels)) {
    intersects[[sprintf("%sx%s", i, j)]] <- gIntersection(Soils[i,],
                                                          Parcels[j,]) 
  }
}

result <- list()
## let's try Parcels, transfer data attributes to new pieces
for (i in 1:nrow(Parcels)) {
  for (j in seq_along(intersects))
   if(gContains(Parcels[i,], intersects[[j]])) {
     result <- c(result, SpatialPolygonsDataFrame(intersects[[j]],     as.data.frame(Parcels[i,]), match.ID = FALSE))

   }
}


## plot
plot(Parcels, xlim = range(c(bbox(Parcels)[1,], bbox(Soils[1,]))),
     ylim = range(c(bbox(Parcels)[2,], bbox(Soils[2,]))))
plot(Soils, add = TRUE)

cols <- colorRampPalette(c("lightblue", "darkblue"))(length(result))
for (i in 1:length(result)) plot(result[[i]], col = cols[i], add = TRUE)
for (i in 1:length(result)) text(coordinates(result[[i]]), label =     as.data.frame(result[[i]])[,"soilType"])

Answer 4

Here's the basic idea (do this as a nested loop over parcels then over soils; I'm not sure it can be vectorized the way that the g* functions are written) :

i <- 2
j <- 2
pieces <- list()
pieces[["int"]] <- gIntersection(Parcels[i,],Soils[j,])
pieces[["diff1"]] <- gDifference(Parcels[i,],Soils[j,])
pieces[["diff2"]] <- gDifference(Soils[j,],Parcels[i,])
plot(Parcels)
plot(Soils,add=TRUE)
plot(pieces[["int"]],col="red",add=TRUE)
plot(pieces[["diff1"]],col="blue",add=TRUE)
plot(pieces[["diff2"]],col="green",add=TRUE)

That should be enough to get you started. The rest is just looping while keeping track of the pieces and merging them all into one big SPDF.

An alternative approach that's more vectorized is:

pieces <- list()
pieces[["int"]] <- gIntersection(Parcels,Soils,byid=TRUE)
pieces[["diff1"]] <- gDifference(Parcels,Soils,byid=TRUE)
pieces[["diff2"]] <- gDifference(Soils,Parcels,byid=TRUE)

This gives you more pieces than actually exist, for some reason. Then you'll have to go back and rbind them and cull the extra pieces that are gEquals.