4

Using the sample data below, how can I generate rasters and spatial points plot with the same colorkey as in the "manually" joined plot shown below?

library(rasterVis)
library(raster)
library(colorRamps)
col=colorRampPalette(matlab.like2(255))

s <- stack(replicate(2, raster(matrix(runif(100), 10))))
xy <- data.frame(coordinates(sampleRandom(s, 10, sp=TRUE)),
                 z1=runif(10), z2=runif(10))

levelplot(s, margin=FALSE, at=seq(0, 1, 0.05),col.regions=col)
x=xy$x;y=xy$y;z=xy$z1

levelplot(z ~ x + y,contour=F, panel = panel.levelplot.points, 
          margin=FALSE,col.regions=col,
          par.settings=list(axis.line=list(lwd=3), strip.border=list(lwd=3)),
         cex=1.4, scales=list(x=list(cex=1.7),y=list(cex=1.7)),xlab=list(label="Longitude",cex=2),
          ylab=list(label="Latitude",cex=2))

sample plot

Thanks to @fdestch I was able to generate the following plot using:

latticeCombineGrid(mget(rep("pp", 24)), layout = c(3, 8))

following my comments on printing multiple plots with the same colorkey.

An issue that remains to be clarified:

1) How can one decide on the order of panels? That is, which row & column to place a particular plot just as in levelplot using index.cond.

enter image description here

code123
  • 2,082
  • 4
  • 30
  • 53

2 Answers2

4

First of all, you should probably make sure that the breaks in the points plot are identical with those defined in the first levelplot.

## raster plot with colorkey disabled
pr <- levelplot(s, margin = FALSE, at = seq(0, 1, 0.05), col.regions = col, 
                colorkey = FALSE, xlab = list("Longitude", col = "transparent"))

## points plot 
pp <- levelplot(z ~ x + y, panel = panel.levelplot.points, cex = 1.4, 
                contour = FALSE, margin = FALSE, col.regions = col, 
                colorkey = list(at = seq(0, 1, .05), width = .6, height = .6),
                xlab = "Longitude", ylab = "Latitude")

Please note the definition of a transparent xlab when creating the raster plot. This little workaround comes in quite handy when using downViewport later on to ensure that the actual plot boundaries of pr and pp overlap (feel free to run grid.rect() right after print(pr, newpage = FALSE) to see what I mean).

The actual plot arrangement can then easily be achieved by using viewports from the grid package.

library(grid)
library(lattice)

## initialize new grid device
grid.newpage()

## add raster plot
vp1 <- viewport(x = 0, y = 0, width = .5, height = 1, 
                just = c("left", "bottom"))

pushViewport(vp1)
print(pr, newpage = FALSE)

## add points plot
downViewport(trellis.vpname("page"))

vp2 <- viewport(x = 1, y = 0, width = .75, height = 1, 
                just = c("left", "bottom"))
pushViewport(vp2)
print(pp, newpage = FALSE)

arranged_plots

fdetsch
  • 5,239
  • 3
  • 30
  • 58
  • very nice solution. Thanks a lot. What if I have 8 rows and 3 columns to display with a single colorkey? Usually I will use `layout=c(3, 8),index.cond=list(c( 1,2,....))` to get my grid. `Viewport` can be tricky for arranging plots. Any work-around for `layout` and `index.cond`? – code123 Apr 07 '16 at 01:25
  • For this particular purpose, there is a function called [`latticeCombineGrid`](https://github.com/environmentalinformatics-marburg/Rsenal/blob/master/R/latticeCombineGrid.R) in our **Rsenal** package. Simply install the package via `devtools::install_github("environmentalinformatics-marburg/Rsenal")` and run e.g. `latticeCombineGrid(mget(rep("pp", 24)), layout = c(3, 8))`. Note that the results look best when all plots have identical x and y limits. – fdetsch Apr 07 '16 at 06:10
  • you saved my day. `latticeCombineGrid` looks very promising. I will definitely try it out. Thanks again. – code123 Apr 07 '16 at 14:42
  • For future reference, `latticeCombineGrid()` has been moved to [Orcs](https://rdrr.io/cran/Orcs/man/latticeCombineGrid.html). – fdetsch Aug 13 '20 at 09:23
  • In above plot, how x and y tick labels can be printed as Geographic coordinates (e.g. 101W, 102W, 15N, etc.)? – raghav Sep 10 '20 at 20:53
  • The raster object needs a valid CRS assigned in order to render geographic coordinates correctly. Something like `projection(s) = "+init=epsg:4326"` should do the trick. For the points data, this basically also applies, but the coding involved becomes a little more extensive. I suggest you have a look at `spplot()` or **sf**'s `plot()` method to figure out this one. – fdetsch Sep 11 '20 at 06:42
1

Here is my solution using latticeExtra::c.trellis:

library(raster)
library(rasterVis)

s <- stack(replicate(2, raster(matrix(runif(100), 10))))
xy <- data.frame(coordinates(sampleRandom(s, 10, sp=TRUE)),
                 z1=runif(10), z2=runif(10))

## Define theme and breaks
myTheme <- BTCTheme()
my.at <- seq(0, 1, 0.05)
  • Plot the Raster* object, using rasterVis::levelplot:

    p1 <- levelplot(s, margin=FALSE,
                        at = my.at,
                        par.settings = myTheme)
    
  • Plot the points, using lattice::levelplot:

    p2 <- levelplot(z1 ~ x + y, data = xy,
                        at = my.at, 
                        panel = panel.levelplot.points,
                        par.settings = myTheme)
    
  • Join them with latticeExtra::c.trellis:

    p3 <- c(p1, p2, layout = c(3, 1))
    
  • Unfortunately, c.trellis does not assign the strip labels correctly, so you have to define them directly:

    update(p3,
        strip = strip.custom(factor.levels = c(names(s), "Points")))
    

raster + points

Oscar Perpiñán
  • 4,491
  • 17
  • 28
  • this is even more straightforward and easy to implement. Thanks. The last line of code can be used to decide where each plot is placed in the panel. If I know the names of `s` then I can easily decide on the order in which to display them. – code123 Apr 09 '16 at 18:13