Is it possible to conditionally subset (and summarize) a raster using another raster as an index?
I have a set of rasters that show monthly data across a number of years. For each point, I need to get the average of just part of the year (e.g., months 2-4). This is very simple using terra::tapp:
library(terra)
#sample data - only five layers for simplicity
s <- rast(nrow=4, ncol=3, nlyrs=5)
n <- ncell(s)
values(s[[1]]) <- sample(1:100, n, replace=T)
values(s[[2]]) <- sample(1:10, n, replace=T)
values(s[[3]]) <- sample(1:10, n, replace=T)
values(s[[4]]) <- sample(1:10, n, replace=T)
values(s[[5]]) <- sample(1:100, n, replace=T)
#take average across layers; separate average of first month from average of months 2-4
m <- tapp(s, index=c(1,2,2,2,1), fun=mean)
#second layer of output gives months 2-4
plot(m[[2]])
The complication is that the exact layers I need vary depending on the pixel. For example, for some pixels, I need to take layers 2-4, but for others I need layers 2-5.
I have a separate index raster (same size/resolution) that gives me, for each pixel, which layer I need to start with, and another that tells me which one I need to end with. Is there any way to conditionally apply the tapp index depending on the cell? Essentially, to create the index using a separate function?
So for the example below, I'd like to have half the cells use idx = c(1,2,2,2,2), while the other half would use idx = c(1,2,2,2,1).
#modify layer 5 to include two different cases of data
values(s[[5]]) <- c(sample(1:10, n/2, replace=T),
sample(1:100, n/2, replace=T))
#index of starting/ending values for the summary statistics
idx_start <- idx_end <- rast(nrow=4, ncol=3)
values(idx_start) <- rep(2, n)
values(idx_end) <- c(rep(4, n/2), rep(5, n/2))
# index for the first raster cell - I'd like this to be responsive across all raster cells
startval <- as.numeric(idx_start[1])
endval <- as.numeric(idx_end[1])
idx <- rep(1, dim(s)[3])
idx[startval:endval] <- 2
idx
#so I could do something like this
m2 <- tapp(s, idx, fun=mean)
You can use terra::rapp (the r stands for "range")
rapp(s, idx_start, idx_end, mean)
#class : SpatRaster
#dimensions : 4, 3, 1 (nrow, ncol, nlyr)
#resolution : 120, 45 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
#source(s) : memory
#name : lyr1
#min value : 3.666667
#max value : 22.250000