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)