2

I am trying to plot 2 rasters with different values with the same scale for comparison. I have been trying to work out the answer using this and this.

I am running into two issues:

  1. I would like to be able to have the legend labels in scientific notation (there are too many decimal places to be of any use)
  2. I would like there to be one scale for both plots, at the bottom.

I am familiar with using legend(), but am not sure how to set the position/direction of the scale and the scientific notation using plot. I have generated some sample data to show the problem I am running into:

library(raster)
set.seed(10)
x = runif(2000000, -0.0005, .9875)
y = runif(2000000, -0.0008, .99)
xmat = matrix(x, nrow = 500)
ymat = matrix(x, nrow = 500)
xras = raster(xmat)
yras = raster(ymat)
min_ = min(minValue(xras), minValue(yras))
max_ = max(maxValue(xras), maxValue(yras))

To plot the data I am using:

par(mfrow = c(2,1))
plot( xras, col=rev( rainbow( 99, start=0,end=1 ) ), breaks=seq(min_,max_,length.out=100), legend = FALSE,  axes = FALSE )
plot( yras, col=rev( rainbow( 99, start=0,end=1 ) ), breaks=seq(min_,max_,length.out=100), legend = FALSE, axes = FALSE )
r.range = c(min_, max_)

And to generate the legend I am using:

plot(xras, legend.only=TRUE, col=rev( rainbow( 99, start=0,end=1 ) ),  legend.width=1, legend.shrink=0.75, axis.args=list(at=seq(r.range[1], r.range[2], abs(r.range[1]-r.range[2])/4), labels=seq(r.range[1], r.range[2], abs(r.range[1]-r.range[2])/4), cex.axis=0.6), legend.args=list(text='Precipitation (m)', side=1, font=2, line=2.5, cex=0.8))

I am open to any method/packages to solve this problem. I have seen that image() is a popular solution for plotting rasters, but would prefer to be able to keep the projection associated with the raster.

To reiterate what I am hoping for: two raster plots displayed (that have different, but similar values (see example)), with a common legend underneath, using scientific notation.

Community
  • 1
  • 1
Sarah
  • 731
  • 7
  • 15
  • You could give a try to the `rasterVis` package. –  Feb 05 '15 at 04:04
  • thanks @Pascal, I've been looking at rasterVis. Trying to figure out how to plot two rasters on one plot (par(mfrow=c(2,1)) seems to be overwritten, and unsure of how to remove the additional plots of the axes. – Sarah Feb 05 '15 at 04:08

2 Answers2

1

This is a try. It will need more refinement:

library(rasterVis)
set.seed(10)
x = runif(2000000, -0.0005, .9875)
y = runif(2000000, -0.0008, .99)
xmat = matrix(x, nrow = 500)
ymat = matrix(x, nrow = 500)
xras = raster(xmat)
yras = raster(ymat)
min_ = min(minValue(xras), minValue(yras))
max_ = max(maxValue(xras), maxValue(yras))
r.range = c(min_, max_)

levelplot(stack(xras, yras), col.regions = rev(rainbow(99, start=0, end=1)), colorkey = list(space = "bottom"))

enter image description here

  • 1
    This worked - the rasters automatically stacked on top of each other (which I didn't expect) - since they are the globe (longer length then height). Thanks! – Sarah Feb 05 '15 at 04:16
  • Maybe I cannot see this clearly but in the legend there seems to be two reds (one above 1.0 and one below 0.0). How do you fix this? – 89_Simple Apr 03 '18 at 13:54
1

Although @Pascal was able to answer the question, I found that I had more flexibility working off my original attempt. Although I said scientific notation, I was incorrect, and just wanted truncated decimal places, this was completed by adding a round function:

axis.args=list(at=seq(r.range[1], r.range[2], abs(r.range[1]-r.range[2])/4),                      labels=round(seq(r.range[1], r.range[2], abs(r.range[1]-.range[2])/4),2),  cex.axis=1),

I was able to put the legend underneath (and modified the inner and outer margins), and adjust settings on it using legend.shrink, legend.width, and the legend.args. In the end, the script was:

library(raster)

set.seed(10)
x = runif(2000000, -0.0005, .9875)
y = runif(2000000, -0.0008, .99)
xmat = matrix(x, nrow = 500)
ymat = matrix(x, nrow = 500)
xras = raster(xmat)
yras = raster(ymat)
min_ = min(minValue(xras), minValue(yras))
max_ = max(maxValue(xras), maxValue(yras))
op <- par(mfrow = c(2,1),
  oma = c(2,4,0,0) + 0.1,
  mar = c(5,0,1,1) + 0.1)
plot( xras, col=rev( rainbow( 99, start=0,end=1 ) ), breaks=seq(min_,max_,length.out=100), legend = FALSE,  axes = FALSE )
plot( yras, col=rev( rainbow( 99, start=0,end=1 ) ), breaks=seq(min_,max_,length.out=100), legend = FALSE, axes = FALSE )
r.range = c(min_, max_) 
plot(xras, col=rev( rainbow( 99, start=0,end=1 ) ), horizontal=TRUE, `breaks=seq(min_,max_,length.out=100), legend.only=TRUE, legend.shrink = 1, legend.width = 3, axis.args = list(at=seq(r.range[1], r.range[2], abs(r.range[1]-r.range[2])/4), labels = round(seq(r.range[1], r.range[2], abs(r.range[1]-r.range[2])/4),2),cex.axis=1), legend.args= list(text='Precipitation (XXX)', side=1, font=3, line = 2, cex=1))`

And looked like this with the sample data: enter image description here

And this with my real data (I changed to axes = TRUE in the plot() commands): enter image description here

Sarah
  • 731
  • 7
  • 15