Overwrite values of raster using highest pixel value across a set of rasters

Question

The following describes a scaled-down example of what I am trying to do.

I have 3 rasters that share the same extent and resolution. The first raster consists of '0' and '1' pixels, the second raster consists of '0' and '2' pixels, and the third raster consists of '0' and '3' pixels. I want to overlay them all together, such that at each pixel, the highest value across the three rasters is preserved. For example, if a pixel has values '1', '2' and '0' for the first, second and third rasters respectively, then the final post-overlay output would have a value '2' for that pixel.

I was initially trying out the raster::overlay and raster::update functions, but did not manage to get them to work. However, after some experimentation, I managed to succeed using the stack(), stackApply() and calc() functions. Although I have managed to get it to work, I am wondering if there are any better alternatives, as eventually, I will be performing this process for 20 rasters that are much larger in extent.

Below is my attempt at a reproducible code as described earlier.

library(raster)

# Create matrix of values to assign to rasters
set.seed(1)
val1 = round(matrix(runif(5*5,0,1),5,5)) # matrix of '0' and '1's
set.seed(2)
val2 = round(matrix(runif(5*5,0,1),5,5))*2 # matrix of '0' and '2's
set.seed(3)
val3 = round(matrix(runif(5*5,0,1),5,5))*3 # matrix of '0' and '3's

# View matrices
val1
val2
val3

# Create random raster and assign matrix values
raster1 = raster(
  nrows = 5,
  ncols = 5,
  xmn = 0,
  xmx = 5,
  ymn = 0,
  ymx = 5,
  vals = val1
) 

raster2 = raster(
  nrows = 5,
  ncols = 5,
  xmn = 0,
  xmx = 5,
  ymn = 0,
  ymx = 5,
  vals = val2
) 

raster3 = raster(
  nrows = 5,
  ncols = 5,
  xmn = 0,
  xmx = 5,
  ymn = 0,
  ymx = 5,
  vals = val3
) 

# View rasters  
plot(raster1)
plot(raster2)
plot(raster3)

# Stack rasters
stack = stack(raster1, raster2, raster3)
plot(stack)

# Create single raster using the maximum values at each location
maximum = stackApply(stack, indices = rep(1, nlayers(stack)),fun = max)
maximum = calc(stack,function(x){max(x)})
plot(maximum)

mrhellmann mrhellmann · Accepted Answer · 2020-10-10T15:59:04

max should work on several different rasters at the same time.

max_raster <- max(raster1, raster2, raster3)

plot(max_raster)

Looks like it's the same as your raster stacked.

stack2 <- stack(maximum, max_raster)
plot(stack2)

A help page with some functions here: http://finzi.psych.upenn.edu/library/raster/html/raster-package.html

Overwrite values of raster using highest pixel value across a set of rasters

1 Answers