8

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
library(raster)
library(ggplot2)

#Create a random raster layer
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))
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
head(rdf)

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"
plot(ndvi)

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()?

2
  • 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
    Commented Oct 24, 2020 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
    Commented Oct 24, 2020 at 20:42

1 Answer 1

11

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
library(rasterVis)
gplot(ndvi) + 
  geom_tile(aes(fill = value)) +
  facet_wrap(~ variable) +
  scale_fill_gradientn(colours = rev(terrain.colors(225))) +
  coord_equal()
  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))) +
  coord_equal()

enter image description here

3
  • 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
    Commented Oct 26, 2020 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
    Commented Oct 26, 2020 at 14:49
  • Ok thanks. This helps.
    – SA Khan
    Commented Oct 29, 2020 at 0:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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