如何执行两个 SpatialPolygonsDataFrame 对象的矢量叠加？

我有两个 GIS 图层 - 称呼它们为Soils and Parcels-- 存储为SpatialPolygonsDataFrames (SPDFs），我想“覆盖”它们，在此处描述的意义上 http://resources.esri.com/help/9.3/ArcGISengine/java/Gp_ToolRef/geoprocessing/overlay_analysis.htm.

覆盖操作的结果应该是一个新的 SPDF，其中：

1. The SpatialPolygons组件包含由两层相交形成的多边形。 （想象一下通过在高射投影仪上叠加两个聚酯薄膜而形成的所有原子多边形）。

2. The data.frame组件记录了组件的属性Soils and Parcels每个原子多边形所属的多边形。

我的问题：是否有现有的 R 函数可以执行此操作？ （我什至很乐意了解一个只得到SpatialPolygons组件右，计算由两层相交形成的原子多边形。）我觉得rgeos应该有一个至少执行（1）的函数，但似乎没有......

这是一个可能有助于让我更清楚我所追求的内容的图，后面是创建Soils and Parcels层数如图所示。

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

自2014年1月起，raster套装包括一个union()使这变得简单的函数：

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)