用分类变量的水平图覆盖栅格的水平图面板

Question

我想覆盖levelplot() 制作的图像。在光栅的每个面板上叠加 16 个格子对象的列表是一场噩梦。我可以使用网格资源单独覆盖地图，但我更喜欢 levelplot() 的美感。让我为如何实现地图制作一个可重现的代码。

lapply(c("raster", "rasterVis", "RColorBrewer"), require, character.only = TRUE)  
# create a rasterbrick
my.brick <- vector("list", 16)
for(i in 1:16){
  my.brick[[i]] <- raster(xmn = 30, xmx = 42, ymn = -6, ymx = 6)
  values(my.brick[[i]]) <- runif(ncell(my.brick[[i]]), 1, 250)
}
my.brick <- stack(my.brick) # the rasterlayers

cols <- c("white", brewer.pal(9, "Reds")) # my colours
my.at <- seq(0, 250, 25)

p1 <- levelplot(my.brick, col.regions = cols, at = my.at) # plot the raster
coods <- vector("list", length = 16) # object to hold my point data
my.plots <- coods  # An empty list to hold trellis objects

for(i in 1:16){ # A loop to create the spatial points
  coods[[i]] <- data.frame(lat = runif(10, 30, 50), lon = runif(10, 0, 20), 
                           val = runif(10, 1, 250))
  attach(coods[[i]])  ## let the headers become R objects
  coordinates(coods[[i]]) <- ~ lon + lat # convert to a sp object.
  my.plots[[i]] <- levelplot(val ~ lon + lat, col.regions = cols, at = my.at, 
                             panel = panel.levelplot.points, cex = 1.3) + 
    layer(sp.points(coods[[i]], pch = 21, bg = "white", col = "black", 
                    lwd = 2, cex = 1.5)) # Plots. The layer is meant enhance 
                                         # the width of the symbol outline
}

如何将my.plots 叠加到光栅图p1 的列表中？ 在堆栈交换上似乎有一些关于此的答案，但我发现它们还不够。

Answer 1

这个问题可以通过previous question 中发布的解决方案来解决：使用latticeExtra::layer 结合panel.number。此外，我修改了您的代码以简化对象创建过程：

library(raster)
library(rasterVis)

## Create a RasterStack
r <- raster(xmn = 30, xmx = 42, ymn = -6, ymx = 6)
my.brick <- stack(lapply(1:16, FUN = function(i)
    init(r,
         fun = function(x) runif(x, min = 1, max = 250)
         )))

## Create a list of SpatialPoints
prj <- CRS(projection(r))
pts <- lapply(1:16, FUN = function(i)
{
    pt <- SpatialPointsDataFrame(cbind(runif(10, 30, 42),
                                       runif(10, -6, 6)),
                                 data.frame(val = runif(10, 1, 250)),
                                 proj4string = prj)
})                    



## Display the raster and the points
levelplot(my.brick) +
    layer(sp.points(pts[[panel.number()]],
                    pch = 21, bg = "white", col = "black", 
                    lwd = 2, cex = 1.5))