2
votes

I've looked through all SO but haven't been able to find an answer to this specific question:

I have a rasterStack with several layers whose values span over a quite large range. Knowing the values of each layer and the chosen colour scale I believe I have managed to plot all the rasters with the same colour scale, but I face three problems now:

  1. I can't be sure that the values are being plotted with the same colour scale, although it seems like it
  2. I can't manage to plot the same scale and scale labels for all layers.
  3. One of my cells in ras3, ras3[3,9], doesn't get coloured, but it is not NA!

Example script follows:

library(raster)
# generate example rasters
set.seed(123)
ras1 <- raster(ncol = 10, nrow= 10)
values(ras1) <- runif(100, 1, 10)

ras2 <- raster(ncol = 10, nrow = 10)
values(ras2) <- runif(100, 5, 50)

ras3 <- raster(ncol = 10, nrow = 10)
values(ras3) <- runif(100, 10, 100)

# stack them 
rasStack <- stack(ras1, ras2, ras3)
# plot normally to check the values
plot(rasStack)

# obtain max and min values 
maxv <- max(maxValue(rasStack))+1
minv <- min(minValue(rasStack))

# set the breaks between min and max values
brks <- seq(minv,maxv,by=0.1)
nbrks <- length(brks)-1
r.range <- c(minv, maxv)

# generate palette
colfunc<-colorRampPalette(c("springgreen", "royalblue", "yellow", "red"))

# plot in a loop with a common legend, using legend.only = T    
for(i in seq_len(nlayers(rasStack))){
  tmp <- rasStack[[i]]
  plot(tmp, breaks=brks,col=colfunc(nbrks), legend = F, zlim=c(minv,maxv), 
  main = names(tmp))
  plot(tmp, legend.only=TRUE, col=colfunc(nbrks),
       legend.width=1, legend.shrink=0.75,
       axis.args=list(at=seq(minv, maxv, 5),
                  labels=round(seq(r.range[1], r.range[2], 5), 2), 
                  cex.axis=0.6),
   legend.args=list(text='value', side=4, font=2, line=2.5, cex=0.8))}

# check that the blank cell has a valid value
ras3[3, 9]
> 99.01704 

Any help will be greatly appreciated!


EDIT: as per ycw's answer I've edited the code and now problem no. 3 has disappeared!

2

2 Answers

1
votes

The last break number should be larger than the maximum value of your data (maxv) so that the cell with the maximum can be colored based on the last color category. Otherwise, the cell will be blank.

I modified your code by changing maxv <- max(maxValue(rasStack)) + 1 but did not change other parts. The result looks good.

2
votes

I just fixed this problem so I'll post the solution in case someone else stumbles with this.

I might be a bit of a workaround, and it certainly is not very elegant, but it works:

First of all we add up all three raster layers in a new one

rasTot <- ras1 + ras2 + ras3

Now we run the loop from the previous code but in the plot with the legend.onlycall we use this total raster.

for(i in seq_len(nlayers(rasStack))){
  tmp <- rasStack[[i]]
  plot(tmp, breaks=brks,col=colfunc(nbrks), legend = F, zlim=c(minv,maxv), 
  main = names(tmp))
  plot(rasTot, legend.only=TRUE, col=colfunc(nbrks),
       legend.width=1, legend.shrink=0.75,
       legend.args=list(text='value', side=4, font=2, line=2.5, cex=0.8))
}

I also edited out some of the legend label specifications, as the defaults are OK.