在由levelplot生成的栅格地图中添加XY点

Question

我有使用R语言的raster包生成的raster地图。可以使用rasterVis包中的levelplot函数来可视化这些栅格地图。

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

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

应该添加每个地图的点，如果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生成的地图不对齐。

- code123

Answer 1

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)

请注意，在layer函数中，columns参数（还有rows参数）指定要将图层添加到哪个面板。

- jbaums

---

完美！你刚刚为我节省了世界的时间。 - code123

现在应该可以正常工作了@JerryN。感谢您的提醒。 - jbaums

Answer 2

所以我有一个包含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)

- code123