Reputation: 17
Is there a way to add a different polygon to a raster stack plot built using levelplot? For example, I plotted a raster stack comprised of 6 rasters. I would like to add a different polygon to each panel. However, when I add a layer, that layer is repeated across all six panels & my attempts to add addition polygons are result in those also being added to each panel.
#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
Upvotes: 2
Views: 2418
Reputation: 6661
You can use ggplot2
with rasterVis
, this will give you more flexibility. The facet_wrap
associated with the correct parameter will allow to separate polygons between raster layers views. To do so, you will have to modify the id
column of your fortified spatial polygon to match raster variable
column.
Here is a reproducible example:
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)
Upvotes: 1
Reputation: 3694
Here's a silly example using the packets
argument in latticeExtra::layer
. If you wanted to, you could just add layer after layer like this.
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)
Upvotes: 3