How to downsample/aggregate R rasters with custom function to produce multiple output rasters?

Question

I have a rather complex function that I'm using to compute aggregates of raster-cells (downsampling). In this function I'm distuingishing between multiple cases, where the case determines which summary function is used. As a simplified example (let's assume, 42, 77 and 123 are dynamically computed values dependent on vec):

my_fun <- function(vec, na.rm=TRUE, VERBOSE = FALSE){
  if(max(vec) > 1){
    print("case 1")
    return(42)
  } else if(max(vec) > 2){
    print("case 2")
    return(123)
  } else {
    print("case 3")
    return(77)
  }
}

vhm_2 = raster::aggregate(vhm, 20, fun=my_fun, expand=TRUE, na.rm=FALSE)

This would result in a raster with values 42, 77 and 123. In addition to this I'd like a second raster, indicating which case was used, i.e. a raster saying containing numbers 1-3 for each cell.

Of course I can run it twice, with a second function that just returns the case number (instead of the "computation"), however I'm processing large rasters and the computations (for determining the case, as well as the actual value) are expensive in the first place, so I'd love to avoid this.

The doc explicitly says that fun needs to return a single number, so my hopes that I could simply return both case number AND computed value were shattered. Is there another (computationally efficient) way to aggregate rasters with custom functions where I could get two rasters as output? https://www.rdocumentation.org/packages/raster/versions/3.5-11/topics/aggregate

Answer 1

You can now do this with terra 1.5.15 (currently the development version. You can install it with install.packages('terra', repos='https://rspatial.r-universe.dev')

library(terra)
#terra 1.5.15
r <- rast()
set.seed(1)
values(r) <- sample(100, ncell(r), replace=TRUE)
a <- aggregate(r, 10, range)
a
#class       : SpatRaster 
#dimensions  : 18, 36, 2  (nrow, ncol, nlyr)
#resolution  : 10, 10  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source      : memory 
#names       : lyr.1, lyr.2 
#min values  :     1,    95 
#max values  :     7,   100 

Or with a custom function

a_fun <- function(v){
  m <- max(v)
  if(m < 80){
    c(1,42)
  } else if(m < 95){
    c(2,123)
  } else {
    c(3, 77)
  }
}

a <- aggregate(r, 5, fun=a_fun)
a
#class       : SpatRaster 
#dimensions  : 36, 72, 2  (nrow, ncol, nlyr)
#resolution  : 5, 5  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source      : memory 
#names       : lyr.1, lyr.2 
#min values  :     1,    42 
#max values  :     3,   123 

But in this case you could also opt to do:

a_fun <- function(v){
  m <- max(v)
  if(m < 80){
    1
  } else if(m < 95){
    2
  } else {
    3
  }
}

a <- aggregate(r, 5, fun=a_fun)
b <- classify(a, matrix(c(1,42, 2,123, 3, 77), ncol=2, byrow=TRUE))