[R-sig-Geo] varying polygon layer with panel in spplot

[Matthew Landis]

Dear R-sig-geo followers:

I wonder if anyone can help me with the following problem.  I'd like to 
plot a series of rasters overlaid with polygons (or points), where the 
polygon layer is specific to each panel.  here's an example of my 
nearest attempt using spplot on a Meuse-based dataset (details follow at 
the end). (If someone can do this with ggplot, that's fine with me 
too).  'sgdf' is a SpatialGridDataFrame and pts.list is a list where 
each element is a SpatialPointsDataFrame.

spplot(sgdf, c('cd','cu', 'pb', 'zn'),
        names.attr = c('Cadmium', 'Copper', 'Lead', 'Zinc'),
        as.table = TRUE, main = "Testing",
        panel = function(x, y, z, subscripts, ...){
            panel.gridplot(x, y, z, subscripts, ...)
            sp.points(pts.list[[subscripts]])
        }
)

When I send this, I get SpatialGridDataFrame plotted properly, but the 
points are not.  Instead I get the messages (by panel): Error using 
packet 1\nrecursive indexing failed at level 2; Error using packet 2\nno 
such index at level 1; Error using packet 3\nno such index at level 1; 
Error using packet 4\nno such index at level 1.

If I specify the subscript manually, e.g. "sp.points(pts.list[[1]])", it 
works fine (but of course the points are the same in all panels).

------------------
Here's the code I used to produce the data:

library(sp)
library(lattice) # required for trellis.par.set():
trellis.par.set(sp.theme()) # sets color ramp to bpy.colors()

data(meuse)
coordinates(meuse)=~x+y

library(gstat, pos = match(paste("package", "sp", sep=":"), search()) + 1)
data(meuse.grid)
coordinates(meuse.grid) = ~x+y
gridded(meuse.grid) = TRUE

v.cd = variogram(log(cadmium) ~ 1, meuse)
ok.model.cd = fit.variogram(v.cd, vgm(1, 'Exp', 500, 1))
cd.ok <- krige(log(cadmium) ~ 1, meuse, meuse.grid, model = ok.model.cd)

v.cu <- variogram(log(copper) ~ 1, meuse)
ok.model.cu <- fit.variogram(v.cu, vgm(1, 'Exp', 500, 1))
cu.ok <- krige(log(copper) ~ 1, meuse, meuse.grid, model = ok.model.cu)

v.pb <- variogram(log(lead) ~ 1, meuse)
ok.model.pb <- fit.variogram(v.pb, vgm(1, 'Exp', 500, 1))
pb.ok <- krige(log(lead) ~ 1, meuse, meuse.grid, model = ok.model.pb)

v.zn = variogram(log(zinc)~1, meuse)
ok.model.zn = fit.variogram(v.zn, vgm(1, "Exp", 500, 1))
zn.ok = krige(log(zinc)~1, meuse, meuse.grid, model = ok.model.zn)

sgdf = cd.ok

sgdf[['cd']] <- cd.ok[['var1.pred']]
sgdf[['cu']] <- cu.ok[['var1.pred']]
sgdf[['pb']] <- pb.ok[['var1.pred']]
sgdf[['zn']] <- zn.ok[['var1.pred']]

pts.list <- list(meuse[, 'cadmium'], meuse[, 'copper'], meuse[, 'lead'], 
meuse[, 'zinc'])

spplot(sgdf, c('cd','cu', 'pb', 'zn'),
        names.attr = c('Cadmium', 'Copper', 'Lead', 'Zinc'),
        as.table = TRUE, main = "Testing",
        panel = function(x, y, z, subscripts, ...){
            panel.gridplot(x, y, z, subscripts, ...)
            sp.points(pts.list[[subscripts]])
        }
)

-----------------------
And, FWIW, the session info:

 > sessionInfo()
R version 2.12.2 (2011-02-25)
Platform: x86_64-pc-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United 
States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] maptools_0.7-38 foreign_0.8-42  maps_2.1-5      lattice_0.19-17
[5] sp_0.9-78       gstat_0.9-77

loaded via a namespace (and not attached):
[1] grid_2.12.2  tools_2.12.2

-- 
~~~~~~~~~~~~~~~~~~~~~~~~~~
Matthew Landis, Ph.D.
Research Scientist
ISciences, LLC
61 Main St. Suite 200
Burlington VT 05401
802.864.2999
~~~~~~~~~~~~~~~~~~~~~~~~~~