23

I'm creating some maps from raster files using the "raster" package in R. I'd like to create comparison rasters, showing several maps side by side. It's important for this that the colour scales used are the same for all maps, regardless of the values in each map. For example, if map 1 has values from 0-1, and map 2 has values from 0-0.5, cells with a value of 0.5 should have the same colour on both maps.

For example:

  • map 1 has values from 0 to 1
  • map 2 has values from 0 to 0.5
  • the colour goes from red (lowest) to green (highest)

I would like a value of 0.5 to have the same colour in both maps (i.e. yellow, as halfway between red and green). The current behaviour is that it is yellow in map 1, and green in map 2.

I can't find a way to make this work. I can't see any way to set the range of pixel values to use with the plotting function. setMinMax() doesn't help (as 'plot' always calculates the values). Even trying to set the values by hand (e.g. g1@data@max <- 10) doesn't work (these are ignored when plotting).

Finally, making a stack of the maps (which might be expected to plot everything on the same colour scale) doesn't work either - each map still has it's own colour scale.

Any thoughts on how to do this?

EDIT:

The solution I ended up using is:

plot( d, col=rev( rainbow( 99, start=0,end=1 ) ), breaks=seq(min(minValue( d )),max(maxValue(d)),length.out=100) )
mo-seph
  • 6,073
  • 9
  • 34
  • 35
  • the solution you ended up using is actually the best answer. Please put it as an answer and accept!! I almost missed it - I don't look for answer in a question. – Tomas Aug 26 '13 at 08:09
  • 1
    `spplot` displays a `RasterStack` with a global legend for all the layers. The `raster` package defines a method for `spplot`, so you don't need to convert your `RasterStack` to an `Spatial*` object to use it. – Oscar Perpiñán Aug 28 '13 at 22:44

6 Answers6

17

Easy solution now is to use the zlim option. 

plot( d, col=rev( rainbow( 99, start=0,end=1 ) ),zlim=c(0,1) )
glaroc
  • 453
  • 4
  • 11
12

Since the image::raster function specifies that the image::base arguments can be passed (and suggests that image::base is probably used), wouldn't you just specify the same col= and breaks= arguments to all calls to image::raster? You do need to get the breaks and the col arguments "in sync". The number of colors needs to be one less than the number of breaks. The example below is based on the classic volcano data and the second version shows how a range of values can be excluded from an image:

 x <- 10*(1:nrow(volcano))
 y <- 10*(1:ncol(volcano))
 image(x, y, volcano, col = terrain.colors( length(seq(90, 200, by = 5))-1), axes = FALSE, breaks= seq(90, 200, by = 5) )
 axis(1, at = seq(100, 800, by = 100))
 axis(2, at = seq(100, 600, by = 100))
 box()
 title(main = "Maunga Whau Volcano", font.main = 4)



 x <- 10*(1:nrow(volcano))
 y <- 10*(1:ncol(volcano))
 image(x, y, volcano, col = terrain.colors( length(seq(150, 200, by = 5))-1), axes = FALSE, breaks= seq(150, 200, by = 5) )
 axis(1, at = seq(100, 800, by = 100))
 axis(2, at = seq(100, 600, by = 100))
 box()
 title(main = "Maunga Whau Volcano Restricted to elevations above 150", font.main = 4)

A specific example would aid this effort.

IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • That's perfect, thanks so much! I ended up doing plot( d, col=rev( rainbow( 99, start=0,end=1 ) ), breaks=seq(min(minValue( d )),max(maxValue(d)),length.out=100) ) (sorry, can't seem to format in the comments!) – mo-seph Nov 30 '10 at 18:31
  • Edit your question and save your solution for posterity. It might come handy someday. – Roman Luštrik Nov 30 '10 at 18:40
5

Added as an answer in response to @Tomas

The answer I ended up using is:

plot( d, col=rev( rainbow( 99, start=0,end=1 ) ), 
    breaks=seq(min(minValue( d )),max(maxValue(d)),length.out=100) )
mo-seph
  • 6,073
  • 9
  • 34
  • 35
4

There is more work to be done here in 'raster' but here is a hack:

 library(raster)
 r1 <- r2 <- r3 <- raster(ncol=10, nrow=10)
 r1[] <- runif(ncell(r1))
 r2[] <- runif(ncell(r2)) / 2
 r3[] <- runif(ncell(r3)) * 1.5
 r3 <- min(r3, 1)
 s <- stack(r1, r2, r3)


 brk <- c(0, 0.25, 0.5, 0.75, 1)
 par(mfrow=c(1,3))
 plot(r1, breaks=brk, col=rainbow(4), legend=F)
 plot(r1, breaks=brk, col=rainbow(4), legend.only=T, box=F)
 plot(r2, breaks=brk, col=rainbow(4), legend=F)
 plot(r1, breaks=brk, col=rainbow(4), legend.only=T, box=F)
 plot(r3, breaks=brk, col=rainbow(4), legend=F)
 plot(r1, breaks=brk, col=rainbow(4), legend.only=T, box=F)

You can also use the spplot function (sp package)

 s <- stack(r1, r2, r3) 
 sp <- as(s, 'SpatialGridDataFrame')
 spplot(sp)

You can also send the values to ggplot (search the r-sig-geo archives for examples) If your RasterLayer links to a very large file, you might first do, before going to ggplot

r <- sampleRegular(r, size=100000, asRaster=TRUE)

and then perhaps

m <- as.matrix(r)
rrrrrh
  • 41
  • 2
1

It did not work for me. I used this script to split the color scale and select the one more suitable according to my data: 

plot(d, col=rev(heat.colors(8, alpha = 1)), breaks = seq(0, 0.40, by = 0.05))
tuomastik
  • 4,559
  • 5
  • 36
  • 48
Kita Romero-Álvarez
  • 11
  • 1
1

A pretty simple solution that should usually work (e.g. with the "plot" function in the raster package) is to set "z axis" limits (which control the colors and the color legend).

E.g. you can do something like: plot(d, zlim=c(0,1))

where d is a stacked raster object. Or, if you have a bunch of separate rasters d1, d2, d2..., you can just do: plot(d1, zlim=c(0,1)) plot(d2, zlim=c(0,1)) plot(d3, zlim=c(0,1)) ...

Adam Clark
  • 86
  • 1
  • 4