[R-sig-Geo] Getting vector index (i.e. layer number) for each cell in raster brick

[Lyndon Estes]

Dear List,

I was hoping to get some advice on the best way to solve the following
problem given two raster bricks.

Given two bricks:

r.list <- lapply(1:20, function(x) {
 r <- raster(xmn = 25, xmx = 32, ymn = -34, ymx = -27)
 res(r) <- 0.1
 r[] <- runif(ncell(r))
 return(r)
})
rs1 <- stack(lapply(1:20, function(x) r.list[[x]]))
rb1 <- brick(rs1)  # Brick 1

r.list2 <- lapply(1:20, function(x) {
 r <- raster(xmn = 25, xmx = 32, ymn = -34, ymx = -27)
 res(r) <- 0.1
 r[] <- sample(1:7000, ncell(r), replace = T)
return(r)
})
rs2 <- stack(lapply(1:20, function(x) r.list2[[x]]))
rb2 <- brick(rs2)  # Brick 2

I am trying to figure out the best way to find, on a per pixel basis,
which layers of brick 1 contains the 10th percentile value, and then
retrieve the value for those layer from another brick 2.

I have gotten as far as using this function to get the quantile's
index for each pixel:

qfun <- function(x) {
  ind <- unname(which(x == quantile(x, 0.1, na.rm = T, type = 3)))
  ind[sample(length(ind))[1]]  # Get rid of ties
}

ind.r <- calc(rb1, fun = qfun)

The above gives me a raster indicating the layer number I need to
extract from brick 2.  However, I am not entirely sure how to access
the layers from there, or if I am even barking up the right tree.

The other alternative I thought of is this:

rst <- seq(1, nrow(rb1), by = nrow(rb1) / 7)
outvals <- lapply(1:length(rst), function(x) {
  j <- getValuesBlock(rb1, row = rst[x], nrows = 10)
  k <- getValuesBlock(rb2, row = rst[x], nrows = 10)
  inds <- unlist(apply(j, 1, function(y) {
   ind <- which(y == sort(y)[2])
   ind2 <- ind[sample(length(ind))[1]]  # Get rid of ties
   return(ind2)
  }))
  j2 <- j[, 1:ncol(j)][cbind(1:nrow(j), inds)]  # Select the correct
column/layer from j
  k2 <- k[, 1:ncol(k)][cbind(1:nrow(k), inds)]  # And from k
  list(j2, k2)
})

r1.q10 <- r.list[[1]] * 0
r2.q10 <- r.list[[1]] * 0

r1.q10[] <- unlist(lapply(1:7, function(x) outvals[[x]][[1]]))  # 10th
percentile values from brick 1
r2.q10[] <- unlist(lapply(1:7, function(x) outvals[[x]][[2]]))  #
Values from same layers as brick 1 10th %iles in brick 2

plot(r1.q10)
plot(r2.q10)

Although this *appears* to work, it seems to be unnecessarily
complicated for something that I am sure is done in a only a few
lines.

I will greatly appreciate any pointers for how to go about solving
this problem.

Many thanks, Lyndon