Overplotting two SpatialPolygonsDataFrames with spplot
Asked Answered
P

1

9

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
Pernas answered 13/10, 2012 at 13:33 Comment(2)
I'd just use base graphics : plot(foo, col=...) and then add=TRUE works and you can control the colours better. Default spplot colours are not good for polygons (white-on-white is never a good idea).Nestor
@Nestor I've been using col.regions=gray(...). I guess the point is taken and I'll move to base graphics or ggplot and reserve lattice for when I need exploratory trellised visualizations. Just means lots more cutting and legending, which I try to avoid where possible.Pernas
A
12

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

enter image description here

Athwartships answered 13/10, 2012 at 15:39 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.