Overplotting two SpatialPolygonsDataFrames with spplot

Question

I have a bunch of data which I've plotted at the county level, without borders. I'd like to add in state boundaries. I have a state shapefile (polygons), but spplot doesn't seem to have any way to add on top of the previous map. Is there any way to do this without rewriting the panel function to take two SPDFs (which seems pretty specialized for what is likely a problem other people have)?

Here's a reproducible example:

library(sp)
Srs1 = Polygons(list(Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))), "s1")
Srs2 = Polygons(list(Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))), "s2")

county <- SpatialPolygonsDataFrame( SpatialPolygons(list(Srs1,Srs2)), 
                                  data.frame( z=1:2, row.names=c("s1","s2") ) )

SrsA <- Polygons(list(Polygon(cbind(c(3,5,5,1,3),c(3,4,6,5,3)))),"sA")
state <- SpatialPolygonsDataFrame( SpatialPolygons(list(SrsA)),
                                  data.frame( z=1,row.names="sA" ))

spplot( county, zcol="z",col=NA )
spplot( state, add=TRUE ) # Note the add=TRUE does nothing here, but that's the spirit of what I want to accomplish

Answer 1

To overplot using the spplot function, you can use the sp.layout argument. For example, create lists of the appropriate layout items, such as

spCounty <- list("sp.polygons", county, col = NA)
spState <- list("sp.polygons", state)

Then plot, passing the above list items as a list to the sp.layout argument:

# spplot(county, zcol = "z", col = NA, sp.layout = list(spCounty, spState))
# actually, you only need to pass the second layout item to sp.layout
spplot(county, zcol = "z", col = NA, sp.layout = spState)

The x and y limits might not be correct if the two spatial data.frames don't overlap completely. You can correct this if necessary by extracting the appropriate limits from bbox(obj)

For example,

theMin <- pmin(bbox(county)[,1], bbox(state)[,1])
theMax <- pmax(bbox(county)[,2], bbox(state)[,2])

spplot(county, zcol = "z", col = NA, sp.layout = spState,
  ylim = c(theMin[2], theMax[2]), xlim = c(theMin[1], theMax[1]))