title: How to perform a vector overlay of two SpatialPolygonsDataFrame objects?
body: <p>I have two GIS layers -- call them <code>Soils</code> and <code>Parcels</code> -- stored as <code>SpatialPolygonsDataFrame</code>s (<code>SPDF</code>s), and I would like to "overlay" them, <a href="http://resources.esri.com/help/9.3/ArcGISengine/java/Gp_ToolRef/geoprocessing/overlay_analysis.htm" rel="noreferrer">in the sense described here</a>.</p>
  
  <p>The result of the overlay operation should be a new SPDF in which:</p>
  
  <ol>
  <li><p>The <code>SpatialPolygons</code> 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).</p></li>
  <li><p>The <code>data.frame</code> component records the attributes of the <code>Soils</code> and <code>Parcels</code> polygons within which each atomic polygon falls.</p></li>
  </ol>
  
  <p><strong>My question(s):</strong> Is there an existing R function that does this? (I would even by happy to learn of a function that just gets the <code>SpatialPolygons</code> component right, calculating the atomic polygons formed from the intersection of the two layers.) I feel like <strong>rgeos</strong> should have a function that does (1) at least, but it seems not to... </p>
  
  <p>Here is a figure that may help make it clearer what I'm after, followed by code that creates the <code>Soils</code> and <code>Parcels</code> layers shown in the figure.</p>
  
  <p><img src="https://i.stack.imgur.com/Tyn8D.png" alt="enter image description here"></p>
  
  <pre><code>library(rgeos)
  
  ## Just a utility function to construct the example layers.
  flattenSpatialPolygons &lt;- function(SP) {
      nm &lt;- deparse(substitute(SP))
      AA &lt;- 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 &lt;-
  local({
      A &lt;- 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 &lt;- flattenSpatialPolygons(A)
      SpatialPolygonsDataFrame(AA,
             data.frame(soilType = paste0("Soil_", LETTERS[seq_along(AA)]),
                        row.names = getSpPPolygonsIDSlots(AA)))
  })
  
  ## Example Parcels layer
  Parcels &lt;-
  local({
      B &lt;- 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 &lt;- flattenSpatialPolygons(B)
      SpatialPolygonsDataFrame(BB,
             data.frame(soilType = paste0("Parcel_", seq_along(BB)),
                        row.names = getSpPPolygonsIDSlots(BB)))
  })
  </code></pre>