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