Reclassifying RasterStack by function and overlaying in R

Question

I've got a RasterStack composed of a four rasters. See below:

# Reproducible examples
kiwi <- raster(xmn = -180, xmx = 180, ymn = -56, ymx = 84, res = 1)
kiwi[] <- sample(2:35, ncell(kiwi), replace = T)

banana <- raster(xmn = -180, xmx = 179.9999, ymn = -56.00081, ymx = 83.99917, 
res = 1)
banana[] <- sample(2:100, ncell(banana), replace = T)

apple <- raster(xmn = -180, xmx = 180, ymn = -56, ymx = 84, res = 1)
apple[] <- sample(4:100, ncell(apple), replace = T)

mango <- raster(xmn = -180, xmx = 180, ymn = -56, ymx = 84, res = 1)
mango[] <- sample(1:90, ncell(mango), replace = T)

# Raster stack
stacked<-stack(kiwi, banana, apple, mango)

I want to reclassify the RasterStack into four categories (tropical, mediterranean, boreal and savannah) based on some conditions (see below), and then combine into one raster and plot all four categories on a single map.

# Example conditions 
 kiwi < 10 & banana > 20 & mango > 0.5 == "tropical" 
 mango > 15 & apple < 10 == "boreal" 

I think I need to create a function and apply it to overlay but I don't know how to write this function. Does anyone know how to do this?

Answer 1

I would prefer raster::calc for this type of function. This is because you can pass it a raster stack or brick and index the vector based on the layers in the stack. With raster::overlay you use x and y arguments, which makes the problem at hand overly complicated.

Let's replicate your example but, note that I am toning down the dimensions a bit (no need to run a n=50,400 example when n=504 will do). Also note that you cannot use character output (eg., tropical, boreal) for a raster, it must be numerically encoded (eg., tropical=1, boreal=2).

library(raster)

kiwi <- raster(xmn = -180, xmx = 180, ymn = -56, ymx = 84, res = 10)
  kiwi[] <- sample(2:35, ncell(kiwi), replace = TRUE)
banana <- kiwi
  banana[] <- sample(2:100, ncell(banana), replace = TRUE)
apple <- kiwi
  apple[] <- sample(4:100, ncell(apple), replace = TRUE)
mango <- kiwi
  mango[] <- sample(1:90, ncell(mango), replace = TRUE)

Here we create a raster stack object.

fruit <- stack(kiwi, banana, apple, mango)   

Now, we write a function that simply indexes the band in the stack (eg., kiwi = 1, banana = 2, apple = 3, mango = 4). Now you can just think of it in terms of which position in a vector do you need to index (eg., for x=c(6,81,20,26), x[1] would be kiwi and a value of 6). With a function in hand, using a nested ifelse function to reclassify the data, we can pass it to raster::calc.

system.class <- function(x, rm.na, ...) {
  return( ifelse(x[1] < 10 & x[2] > 20 & x[4] < 15, 1,
            ifelse(x[4] > 15 & x[3] < 10, 2, 0)) )
}

( fruit.class <- calc(fruit, system.class) )
  plot(fruit.class)

You can visualize what is happening here a bit better by coercing to an sp object. The @data slot is a dataframe containing the raster data by row. That is to say that each row in the dataframe represents a single pixel across all the layers. Using this format we can use the apply function to pass our function to the data. The result will be a vector of the results. This would represent a single band raster.

f <- as(fruit, "SpatialPixelsDataFrame")
  head(f@data)

( f.class <- apply(f@data, MARGIN=1, FUN=system.class) ) 
  f@data$class <- f.class 

plot(raster(f,layer=5))