My aim is to plot a rasterstack as a single image using ggplot2. First, I apply ggplot on a single layer raster using the following code:

#Load packages

#Create a random raster layer
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))
crs(r) <- "+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
plot(r) #Produce simple plot

rdf <- as.data.frame(r, xy=TRUE) #Convert raster to data.frame
names(rdf)[3] <- 'magnitude' #Name value column

Now I apply gglpot()

ggplot(data = rdf)+
  geom_raster(mapping=aes(x=x, y=y, fill=magnitude))+
  scale_fill_gradientn(colours= rev(terrain.colors(10)), name='Magnitude')

It successfully prints. Now I create a raster stack with 10 layers and call this ndvi.

rr <- lapply(1:10, function(i) setValues(r, runif(ncell(r), min = -10000, max = 10000)))
ndvi <- stack(rr)
crs(ndvi) <- "+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

How could I plot this stack using ggplot2? Do I create a data.frame but then how do I supply fill argument in aes() and facet_wrap()?

  • 1
    How would you imagine that you could use all 10 layers in a single graph? Do you want a set of 10 smaller graphs, or somehow combine all 10 into an image? – Simbamangu Oct 24 '20 at 20:37
  • I would prefer 10 graphs into one image using facet_warp() etc. But it would be interesting to learn how to to produce a set of 10 smaller graphs. – SA Khan Oct 24 '20 at 20:42

You can use scale_fill_gradientn to specify the fill, and facet_wrap should pick the colours from there:

  1. If you have rasterVis package installed, you can use gplot function to create the ggplot object
gplot(ndvi) + 
  geom_tile(aes(fill = value)) +
  facet_wrap(~ variable) +
  scale_fill_gradientn(colours = rev(terrain.colors(225))) +
  1. Otherwise, you can construct the data yourself, which is basely what gplot function does
coords <- xyFromCell(ndvi, seq_len(ncell(ndvi)))
ndvi <- stack(as.data.frame(getValues(ndvi)))
names(ndvi) <- c('value', 'variable')

ndvi <- cbind(coords, ndvi)
ggplot(ndvi) + 
  geom_tile(aes(x, y, fill = value)) +
  facet_wrap(~ variable) +
  scale_fill_gradientn(colours = rev(terrain.colors(225))) +

enter image description here

  • Thanks for the help. This is what I needed. One more thing, how can I plot it layer by layer where each layer is a separate plot? – SA Khan Oct 26 '20 at 14:04
  • 1
    You could apply a function (that plots a single layer) over the layer names. So, each time you use the layer name to extract the corresponding layer and create a ggplot object. This means you will end up with a list of ggplot objects. – Liman Oct 26 '20 at 14:49
  • Ok thanks. This helps. – SA Khan Oct 29 '20 at 0:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.