使用R中的levelplot将不同的多边形添加到栅格堆积图中

Question

有没有一种方法可以将不同的多边形添加到使用levelplot构建的光栅堆叠图中？例如，我绘制了一个由6个光栅组成的光栅堆栈。我想在每个面板上添加一个不同的多边形。但是，当我添加一个层时，该层在所有六个面板上重复&我尝试添加附加多边形时，这些多边形也会被添加到每个面板中。

#plot raster stack
levelplot(rStack, margin=F, xlab="", ylab="", 
par.settings=list(strip.background=list(col="lightgray")), 
names.attr=c("18000BP", "15000BP", "12000BP", "6000BP", "2050", "2090"),
      scales = list(draw=F))

levelplot of raster stack without polygons

#plot raster stack & add polygons
levelplot(rStack, margin=F, xlab="", ylab="", 
par.settings=list(strip.background=list(col="lightgray")), 
names.attr=c("18000BP", "15000BP", "12000BP", "6000BP", "2050", "2090"),
      scales = list(draw=F)) + 
layer(sp.polygons(ice1)) + layer(sp.polygons(ice2))

levelplot with two polygons added, which are repeated across all panels

Answer 1

下面是一个在latticeExtra::layer中使用packets参数的愚蠢示例。如果你愿意，你可以像这样一层又一层地添加。

library(lattice)
library(latticeExtra)
library(sp)

Sr1 = Polygon(cbind(c(2, 4, 4, 1, 2), c(2, 3, 5, 4, 2)))
Sr2 = Polygon(cbind(c(5, 4, 2, 5), c(2, 3, 2, 2)))
Sr3 = Polygon(cbind(c(4, 4, 5, 10, 4), c(5, 3, 2, 5, 5)))
Sr4 = Polygon(cbind(c(5, 6, 6, 5, 5), c(4, 4, 3, 3, 4)), hole = TRUE)

Srs1 = Polygons(list(Sr1), "s1")
Srs2 = Polygons(list(Sr2), "s2")
Srs3 = Polygons(list(Sr3, Sr4), "s3/4")
SpP = SpatialPolygons(list(Srs1, Srs2, Srs3), 1:3)

levelplot(rnorm(10, 10) ~ 1:10 + runif(10, 1, 10) | rep(c("A", "B"),each = 5)) +
  layer(sp.polygons(SpP), packets = 1)

​

​

Answer 2

您可以在rasterVis中使用ggplot2，这将为您提供更大的灵活性。与正确参数关联的facet_wrap将允许在栅格层视图之间分隔多边形。为此，必须修改加强的空间多边形的id列，使其与光栅variable列匹配。

下面是一个可重复使用的示例：

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

# Raster object
fn <- system.file("external/test.grd", package="raster")
s <- stack(fn, fn)
names(s) <- c("A", "B")

# Two polygons
pol.A <- data.frame(x = c(179000, 179000, 180000, 180000),
                    y = c(330000, 331000, 331000, 330000))
pol.B <- data.frame(x = c(180000, 180000, 181000, 181000),
                    y = c(331000, 332000, 332000, 331000))
# Transform as Spatial Polygons to fit your example
pol.sp <- SpatialPolygons(list(Polygons(list(Polygon(pol.A)), "A"),
                                 Polygons(list(Polygon(pol.B)), "B")))
# Use Fortify to be able to use SpatialPolygons with ggplot2
# rename id as "variable" to correspond to rasterstack gplot output
pol.sp.fort <- fortify(pol.sp) %>%
  rename(variable = id)

# Plot. Note that gplot is not ggplot and issued from library rasterVis
gplot(s) + 
  geom_tile(aes(fill = value)) +
  geom_polygon(data = pol.sp.fort, aes(long, lat), col = "red") + 
  facet_wrap(~ variable)