将 XY 点添加到 levelplot 生成的栅格地图

Add XY points to raster map generated by levelplot

我有 raster 地图,它们是使用 R 中的 raster 包生成的。这些栅格层可以使用 rasterVis 包的 levelplot 函数可视化:

levelplot(rasterstack, layout=c(1, 2), 
          col.regions=colorRampPalette(c('darkred', 'red3', 'orange2', 'orange', 
                                         'yellow', 'lightskyblue', 'steelblue3', 
                                         'royalblue3', 'darkblue')))

现在,我想将一些由 xy 坐标定义的 z 值添加到 levelplot 地图。包含 z 值的数据框有 4 列。第 1 列和第 2 列包含 x 和 y 坐标,第 3 列包含布局 (1, 1) 中地图 1 的 z 值和布局 (1, 2) 中的第 4 列。

应添加每张地图的点数,如果 z < 0.05,则 pch=2,如果 z > 0.05,则 pch=3

我在网上搜索并找到了 Ripley 的解决方案,但它不适用于我的情况:

levelplot(rcp852, xlab = "", ylab = "",
          panel = function(x, y, subscripts, ...) {
            panel.levelplot(x, y, subscripts, ...)
            panel.xyplot(topo$x,topo$y, cex = 0.5, col = 1)
          }
)

我尝试了许多其他选项,但这些点与通过 levelplot 生成的地图不一致。

layer用起来很方便:

library(raster)
library(rasterVis)
library(sp)

s <- stack(replicate(2, raster(matrix(runif(100), 10))))
xy <- data.frame(coordinates(sampleRandom(s, 10, sp=TRUE)),
                z1=runif(10), z2=runif(10))
coordinates(xy) <- ~x+y

levelplot(s, margin=FALSE, at=seq(0, 1, 0.05)) + 
  layer(sp.points(xy, pch=ifelse(xy$z1 < 0.5, 2, 3), cex=2, col=1), columns=1) +
  layer(sp.points(xy, pch=ifelse(xy$z2 < 0.5, 2, 3), cex=2, col=1), columns=2)

请注意 layercolumns 参数(rows 也存在)指定要将图层添加到哪个面板。

所以我有一个包含 xy 坐标和许多 z 列的数据框。感谢@jbaums,这是将点添加到我的地图的最终答案:

s <- stack(raster1,raster2)
coordinates(SITES...TTEST) <- ~x+y # SpatialPointsDataFrame
levelplot(s, layout=c(1, 2), 
          col.regions=colorRampPalette(c("darkred", "red3", 
                                         "orange", "yellow", "lightskyblue", "royalblue3", 
                                         "darkblue")),
          at=seq(floor(39.15945) ,ceiling(51.85068), length.out=30), 
          par.strip.text=list(cex=0),scales=list(alternating=FALSE))+
  layer(sp.points(SITES...TTEST, pch=ifelse(SITES...TTEST$Precip_DJF6 < 0.05, 2, 3), cex=2, col=1), rows=1) +
  layer(sp.points(SITES...TTEST, pch=ifelse(SITES...TTEST$Precip_DJF6 < 0.05, 2, 3), cex=2, col=1), rows=2)