[GIS] Plotting multiple raster stacks with rasterVis::gplot and ggplot2::facet_wrap

Question

I'm trying to ggplot2::facet_wrap a stack of rasters, and then call geom_raster to add more rasters on top of the previously mapped raster stack.

Here is a very simple example dataset to re-create what I've tried:

    #Load packages
    library(raster)
    library(ggplot2)
    
    # create a raster 
    set.seed(11)
    r <- raster(nrows = 10, ncols = 10, res = 30, xmn = 267195, xmx = 267375, ymn = 4016985, ymx = 4017135)
    r <- setValues(r, runif(ncell(r), min = -10000, max = 10000))
    
    # create a raster stack 
    rr <- lapply(1:4, function(i) setValues(r, runif(ncell(r), min = -10000, max = 10000)))
    r_stack <- stack(rr)
    crs(r_stack) <- "+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
    prj <- crs(r_stack)
    plot(r_stack)
    
    
    # use gplot w/ facet_wrap
    library(rasterVis)
    gplot(r_stack) + 
      geom_raster(aes(fill = value)) +
      facet_wrap(~ variable) +
      scale_fill_gradientn(colours = rev(terrain.colors(225))) +
      coord_equal()
  

   # create another raster stack 
    ss <- lapply(1:4, function(i) setValues(r, runif(ncell(r), min = -1, max = 1)))
    s_stack <- stack(ss)
    crs(s_stack) <- "+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
    prj <- crs(s_stack)

# plot r_stack and call geom_raster to plot s_stack on top
    gplot(r_stack) + 
      geom_raster(aes(fill = value)) +
      facet_wrap(~ variable) +
      scale_fill_gradientn(colours = rev(terrain.colors(225))) +
      coord_equal() + 
      geom_raster(data=s_stack, aes(x=x,y=y))

#Error: `data` must be a data frame, or other object coercible by `fortify()`, not an S4 object with class RasterStack
#Run `rlang::last_error()` to see where the error occurred.
 

I realize that the new raster would essential cover up the previous raster plot, but the raster data that I'm calling in the 2nd geom_raster call will have different data to display, most of which will be 'gaps' or NA.

Answer 1

Two things:

• geom_raster() takes a data frame, not a raster object. You need to convert your stacks before plotting them.
• It is difficult (but not impossible) to set multiple fill aesthetics in the same ggplot - one approach is documented at https://eliocamp.github.io/codigo-r/2018/09/multiple-color-and-fill-scales-with-ggplot2/ but its not behaviour that's built into ggplot's design, and that is deliberate. The preferred option is to bind the above two dataframes together, but you need to have a plan for any surprise overlaps.

that said,

library(tidyr)
library(dplyr)

r_stack_df <- as.data.frame(r_stack, xy = TRUE) %?% 
  tidyr::pivot_longer(cols = !c(x, y), 
                      names_to = 'variable', 
                      values_to = 'value')

s_stack_df <- as.data.frame(s_stack, xy = TRUE) %>%
  tidyr::pivot_longer(cols = !c(x, y), 
                      names_to = 'variable', 
                      values_to = 'value') %>%
  # for contrast:
  dplyr::filter(variable %in% c('layer.1', 'layer.4'))
  
ggplot() + 
  geom_raster(data = r_stack_df, 
              aes(x = x, y = y, fill = value), alpha = 0.2) +
  facet_wrap(~ variable) +
  scale_fill_gradientn(colours = rev(terrain.colors(225))) +
  coord_equal() + 
  theme_minimal() +
  geom_raster(data = s_stack_df, 
              aes(x = x, y = y, fill = value), alpha = 0.2)

or,

x_stack_df <- 
  rbind(dplyr::slice_sample(r_stack_df, prop = 0.4) %>% mutate(hue = 'red'),
        dplyr::slice_sample(s_stack_df, prop = 0.4) %>% mutate(hue = 'blue'))

ggplot() + 
  geom_raster(data = x_stack_df, 
              aes(x = x, y = y, fill = hue), alpha = 0.4) +
  facet_wrap(~ variable) +
  scale_fill_identity() +
  coord_equal() + 
  theme_minimal()

The above was done with ggplot2 v3.3.5, tidyr v1.1.3, dplyr v1.0.7, and R 4.1.0.

You should perhaps consider package tmap, or use a GIS for this kind of visualisation rather than R.