I have a list of multilayer rasters (e.g., stacks) with matching extents and layer names, and I need to average them to create a single stack retaining the appropriate layers. This can easily be achieved by manually removing each raster from the list and simply using terra::mean(stack.1, stack.2, na.rm = T).
Unfortunately, I have to do this repeatedly with numerous files and need an omnibus solution. I am also aware that I can use terra::rast() on a list of rasters to create a single stack, and could take the mean of this object, but this approach averages all layers resulting in only a single layer, rather than a stack of averaged layers. Similarly, lapply() and do.call() are not helpful because these apply functions to each element in the list, not all elements collectively. I would greatly appreciate any suggestions.
Code sample:
# Structure of stacks
stack.1
class : SpatRaster
dimensions : 776, 720, 6 (nrow, ncol, nlyr)
resolution : 0.6, 0.6 (x, y)
extent : 278110.2, 278542.2, 4752927, 4753393 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 16N (EPSG:32616)
source : pepr.2022.100.tif
names : preds~_temp, preds~_temp, preds~_temp, preds~vapor, preds~vapor, preds~vapor
min values : -2.074884, 29.46503, 7.464921, -118.1709, 2325.729, 984.7896
max values : 1.352281, 43.14127, 11.517185, 128.4751, 2929.794, 1289.9479
# Listed
stack.list <- list(stack.1, stack.2, stack.3)
# Desired result, a stack containing the same layers that have been averaged
mean(stack.1, stack.2, stack.3, na.rm =T)
class : SpatRaster
dimensions : 776, 720, 6 (nrow, ncol, nlyr)
resolution : 0.6, 0.6 (x, y)
extent : 278110.2, 278542.2, 4752927, 4753393 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 16N (EPSG:32616)
source(s) : memory
names : preds~_temp, preds~_temp, preds~_temp, preds~vapor, preds~vapor, preds~vapor
min values : -1.928547, 29.14394, 5.955294, -115.4321, 2044.692, 941.3278
max values : 1.215450, 38.19327, 9.235113, 127.4777, 2525.456, 1166.1364
You can create a SpatRasterDataset and then use terra::app
Example data
r <- rast(system.file("ex/logo.tif", package="terra"))
x <- list(r, r*10, r*100)
Solution
s <- sds(x)
app(s, mean)
#class : SpatRaster
#dimensions : 77, 101, 3 (nrow, ncol, nlyr)
#resolution : 1, 1 (x, y)
#extent : 0, 101, 0, 77 (xmin, xmax, ymin, ymax)
#coord. ref. : Cartesian (Meter)
#source(s) : memory
#names : red, green, blue
#min values : 0, 0, 0
#max values : 9435, 9435, 9435
"manual approach" for comparison
minmax(mean(r, r*10, r*100))
# red green blue
#min 0 0 0
#max 9435 9435 9435