0
votes

I want to plot my rasterlayer (b1_mosaic_diff) with gradient colours: Negative values should turn from yellow into red and positive numbers from light blue into dark blue.

I have a code sample which is not working out unfortunately (instead of gradient red i put heat.colors same as topo.colors for gradient blue:

plot(b1_mosaic_diff, col=ifelse(na.omit(as.data.frame(b1_mosaic_diff))<0, heat.colors(n=5), topo.colors(n=5)), main="Difference between mean & median: '+' representing higher mean")

I don't want to use ggplot, as it is easier to me to stick to the normal plot. Is there a possiblity to do this? Thanks in advance!

1
Reproducible example please: stackoverflow.com/a/5963610/1412059Roland

1 Answers

2
votes

This is a little bit fiddly, but I usually generate n (e.g. 1000) colours with each of two complementary colourRampPalettes (e.g. one made of reds/yellows and one of blues). Afterwards, if the colour key does not have even length either side of zero, I subset the colour vector corresponding to the shorter portion of the colour key, according to it's length as a proportion of the longer end of the key (this is easier to explain with the example below). Finally, I concatenate the two vectors of colours and pass it to rasterVis::levelplot's col.regions argument.

Here's an example:

  1. Load required libraries and create a dummy raster:

    library(raster)
    library(RColorBrewer)
    library(rasterVis)
    
    r <- raster(matrix(runif(100, -5, 10), 10))
    
  2. Create a vector of colours. Here I use the RColorBrewer palettes 'YlOrRd' (which I reverse to make it red > orange > yellow, since the colours are associated with raster values in the order of most negative to most positive), and 'Blues'. These palettes only provide 9 colours each, so we pass the palettes to colorRampPalette in order to generate many intervening (interpolated) colours. Here I generate 1000 colours with each palette - this is more than enough for creating a smooth colour key.

    We want the colour key to span the range of values in the raster dataset, which range from -5 to 10. So, we need to subset the vector of colours that will be used for the shorter arm of the colour key (the negative values, i.e. reds). To be exact, we want to drop the first 500 values of the reds vector, since the <0 portion of the key is half the length of the >0 portion. This leaves us with a vector of 500 orange-yellows (we dropped the reddest of the reds), and a vector of 1000 blues.

    cols <- c(colorRampPalette(rev(brewer.pal(9, 'YlOrRd')))(1000)[501:1000], 
              colorRampPalette(brewer.pal(9, 'Blues'))(1000))
    
  3. Plot the raster with levelplot, which is a bit more flexible than the regular plot method for raster objects. We first create a vector of lower and upper bounds (z-limits), for convenience.

    zrng <- range(pretty(c(minValue(r), maxValue(r))))
    
    levelplot(r, margin=FALSE, at=seq(zrng[1], zrng[2], diff(zrng)/1000),
              col.regions=cols, scales=list(draw=FALSE))
    

    enter image description here

  4. Bonus exercise!

    To automate the vector sub-setting, you can use the following function:

    asym_colours <- function(r, pal1, pal2, n=2000) {
      zrng <- range(pretty(c(minValue(r), maxValue(r))))
      brks <- do.call(seq, as.list(c(max(abs(zrng)) * c(-1, 1), length.out=2*n)))
      c(pal1(n)[(sum(brks < zrng[1]) + 1):n], 
        pal2(n)[1:(sum(brks < zrng[2] & brks > 0) + 1)])  
    }
    

    For example:

    r <- raster(matrix(runif(100, -0.3, 0.1), 10))    
    cols <- asym_colours(r, colorRampPalette(rev(brewer.pal(9, 'YlOrRd'))),
                         colorRampPalette(brewer.pal(9, 'Blues')))
    zrng <- range(pretty(c(minValue(r), maxValue(r))))
    levelplot(r, margin=FALSE, at=seq(zrng[1], zrng[2], diff(zrng)/1000),
              col.regions=cols, scales=list(draw=FALSE))
    

    enter image description here

    or, with a different ramp:

    cols <- asym_colours(r, colorRampPalette(rev(brewer.pal(9, 'GnBu'))),
                         colorRampPalette(brewer.pal(9, 'PuRd')))
    zrng <- range(pretty(c(minValue(r), maxValue(r))))
    levelplot(r, margin=FALSE, at=seq(zrng[1], zrng[2], diff(zrng)/1000),
              col.regions=cols, scales=list(draw=FALSE))
    

    enter image description here