Add XY points to raster map generated by levelplot

问题:

I have raster maps which are generated using the raster package in R. These raster layers can be visualized using the rasterVis package's levelplot function:

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

Now, I would like to add some z values defined by xy cordinates to the levelplot map. The dataframe containing z values has 4 columns. Columns 1 & 2 contain x & y coordinates, column 3 contains z values for map 1 in layout(1, 1) and column 4 for layout(1, 2).

The points per map should be added such that if z < 0.05, pch=2 and if z > 0.05, pch=3.

I have searched the web and found a solution by Ripley but it does not work in my case:

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)
          }
)

I tried many other options but the points do not align with the map generated via levelplot.

回答1:

layer is very convenient for this:

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

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

Note that the columns argument to layer (rows also exists) specifies which panel(s) you want to add the layer to.

回答2:

So I had a dataframe with xy cordinates and so many z columns. This is the final answer to get points added to my map thanks to @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)