Making new categorical raster selected from stack having maximum values for each pixels in R

Question

I have a raster stack having different raster layers with band names (e.g. B1, B2, B3, ….…. Bn). I need to prepare a new categorical raster with band name (for each pixels) having maximum pixel values selected from a stack.

I have used the following code:

library(raster) 
files <- list.files("E:/Red Panda/TIF/", pattern='tif$', full.names=TRUE ) 
Raster_stack <- stack(files) 
pixelMaxBand.ras <- whiches.max(Raster_stack) 

But I found this type of error:

Error in .local(x, ...) : you can use only use this function for an object with less than 10 layers.

With “whiches.max” function for layers greater than 10, can you suggest the appropriate code?

Answer 1

You can pass the base function which.max to raster::calc. The vector length should remain constant because raster::calc does not remove NA's unless you pass it a function that supports an na.rm argument. What you are after is a simple as:

library(raster)

r <- do.call(raster::stack, replicate(6, 
       raster::raster(matrix(runif(1000), 100, 100))))

( rmax <- calc(r, which.max) )  

The terra package allows you to apply which.max directly to the raster object. I recommend starting to move your workflows to terra as it is the eventual replacement for raster and is considerably faster with a better function interface.

library(terra)
r <- rast(r)
( rmax <- which.max(r) ) 

Answer 2

You can use which.max directly

Example data

library(raster)
f <- stack(system.file("external/rlogo.grd", package="raster")) 
s <- stack(f/2, f, f, f)
nlayers(s)
#[1] 12
# add some `NA`s
s[[12]] <- setValues(s[[12]], NA)

(r <- which.max(s))   
#class      : RasterLayer 
#dimensions : 77, 101, 7777  (nrow, ncol, ncell)
#resolution : 1, 1  (x, y)
#extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
#crs        : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source     : memory
#names      : layer 
#values     : 1, 6  (min, max)

As Jeffrey Evens pointed out, you can also do:

rc <- calc(s, which.max)

The reason that you can use whiches.max with no more then 9 layers (it could have been 10) is that it also returns ties:

(rr <- whiches.max(s[[1:9]]))
#class      : RasterLayer 
#dimensions : 77, 101, 7777  (nrow, ncol, ncell)
#resolution : 1, 1  (x, y)
#extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
#crs        : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m #+no_defs 
#source     : memory
#names      : layer 
#values     : 47, 123456789  (min, max)

47 means that layers 4 and 7 both have the maximum value.

Also note that which.max always ignores NAs:

which.max(c(NA, 1:5, NA, 3))
#[1] 6

And with terra:

library(terra)    
f <- system.file("ex/logo.tif", package="terra")
x <- rast(c(f, f, f, f))
x[[1:3]] <- x[[1:3]] / 2
which.max(x)
#class       : SpatRaster 
#dimensions  : 77, 101, 1  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source      : memory 
#name        : which.max 
#min value   :         1 
#max value   :         6