Variation of https://gist.github.com/oscarperpinan/5451682 using SpatialPolygons

stamenPoly.R

library(sp)
library(ggmap)
## latticeExtra must be loaded after ggmap because both ggplot2 and
## latticeExtra define a 'layer' function. We need the definition from
## latticeExtra.
library(latticeExtra)
## We only need maptools to get an example
library(maptools)
SIDS <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
                      IDvar="FIPSNO",
                      proj4string=CRS("+proj=longlat +ellps=clrk66"))

## Get the stamen map. Same as with SpatialPoints
bbPoly <- bbox(SIDS)
gmap <- get_map(c(bbPoly), maptype='watercolor', source='stamen', crop=FALSE)
bbMap <- attr(gmap, 'bb')
latCenter <- with(bbMap, ll.lat + ur.lat)/2
lonCenter <- with(bbMap, ll.lon + ur.lon)/2
height <- with(bbMap, ur.lat - ll.lat)
width <- with(bbMap, ur.lon - ll.lon)

## Use latticeExtra::+.trellis and latticeExtra::layer instead of sp.layout
spplot(SIDS['BIR79']) +
                  layer(grid.raster(gmap,
                                    x=lonCenter, y=latCenter,
                                    width=width, height=height,
                                    default.units='native'),
                        under=TRUE)

stamenPolyWithLayerInFunction.R

## If you need to use layer inside a function, you have to evaluate the code inside "layer({})" 
## or pass the arguments as a list in 'data'.
myFoo <- function(){
    spplot(SIDS['BIR79']) + 
        layer({
            bbPoly <- bbox(SIDS)
            gmap <- get_map(c(bbPoly), maptype='watercolor', source='stamen', crop=FALSE)
            bbMap <- attr(gmap, 'bb')
            latCenter <- with(bbMap, ll.lat + ur.lat)/2
            lonCenter <- with(bbMap, ll.lon + ur.lon)/2
            height <- with(bbMap, ur.lat - ll.lat)
            width <- with(bbMap, ur.lon - ll.lon)
            grid.raster(gmap,
                        x=lonCenter, y=latCenter,
                        width=width, height=height,
                        default.units='native')
        }, under=TRUE)
}

stamenPolyWithLayerInFunction2.R

myFoo2 <- function(){
    bbPoly <- bbox(SIDS)
    gmap <- get_map(c(bbPoly), maptype='watercolor', source='stamen', crop=FALSE)
    bbMap <- attr(gmap, 'bb')
    latCenter <- with(bbMap, ll.lat + ur.lat)/2
    lonCenter <- with(bbMap, ll.lon + ur.lon)/2
    height <- with(bbMap, ur.lat - ll.lat)
    width <- with(bbMap, ur.lon - ll.lon)
    
    spplot(SIDS['BIR79']) + 
        layer({
            grid.raster(gmap,
                        x=lonCenter, y=latCenter,
                        width=width, height=height,
                        default.units='native')
        }, under=TRUE,
              data=list(gmap=gmap,
                  lonCenter=lonCenter, latCenter=latCenter,
                  width=width, height=height))
}

stamenPolyWithLayerInFunction3.R

myFoo3 <- function(poly){
    bbPoly <- bbox(poly)
    gmap <- get_map(c(bbPoly), maptype='watercolor',
                    source='stamen', crop=FALSE)
    bbMap <- attr(gmap, 'bb')
    latCenter <- with(bbMap, ll.lat + ur.lat)/2
    lonCenter <- with(bbMap, ll.lon + ur.lon)/2
    height <- with(bbMap, ur.lat - ll.lat)
    width <- with(bbMap, ur.lon - ll.lon)
    
    spplot(poly) + 
        layer({
            grid.raster(gmap,
                        x=lonCenter, y=latCenter,
                        width=width, height=height,
                        default.units='native')
        }, under=TRUE,
              data=list(gmap=gmap,
                  lonCenter=lonCenter, latCenter=latCenter,
                  width=width, height=height))
}

myFoo3(SIDS["BIR79"])

stamenPolyWithLayerInFunction4.R

myFoo4 <- function(poly, map, ...){
    bbMap <- attr(map, 'bb')
    latCenter <- with(bbMap, ll.lat + ur.lat)/2
    lonCenter <- with(bbMap, ll.lon + ur.lon)/2
    height <- with(bbMap, ur.lat - ll.lat)
    width <- with(bbMap, ur.lon - ll.lon)
    
    spplot(poly, ...) + 
        layer({
            grid.raster(map,
                        x=lonCenter, y=latCenter,
                        width=width, height=height,
                        default.units='native')
        }, under=TRUE,
              data=list(map=map,
                  lonCenter=lonCenter, latCenter=latCenter,
                  width=width, height=height))
}

mapWater <- get_map(c(bbox(SIDS)), maptype='watercolor',
               source='stamen', crop=FALSE)

myFoo4(SIDS["BIR79"], mapWater,
       col.regions=terrain.colors(30))

mapTerrain <- get_map(c(bbox(SIDS)), maptype='terrain',
                      source='stamen', crop=FALSE)

myFoo4(SIDS["BIR79"], mapTerrain,
       col.regions=terrain.colors(30))

Hello Oscar,
very nice code, especially the latticeExtra section was very helpful to me.
But may I ask a question...
If I put that stuff inside a function, even if I use print(spplot...), it does not work. It says that the object gmap is not found. Do you have any idea what is the problem and how to solve it?
Thank you and greetings!

I suspect your problem is related to the latticeExtra::layer mechanism. From its help page:

“Note that the evaluation used in ‘layer’ is non-standard, and can be confusing at first: you typically refer to variables as if inside the panel function (‘x’, ‘y’, etc); you can usually refer to objects which exist in the global environment (workspace), but it is safer to pass them in by name in the ‘data’ argument to ‘layer’. (And this should not to be confused with the ‘data’ argument to the original ‘xyplot’.)”

Watercolor stamen maps are now in jpeg format. Version 2.4 of ggmap fixes this: dkahle/ggmap@d657e1c.