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

[Lyndon Estes]

Hi Robert,

Very many thanks, and I see the two hours I spent on that last night
would have been more profitably applied in the help files.

Cheers, Lyndon

On Fri, Dec 2, 2011 at 12:10 AM, Robert J. Hijmans <r.hijmans at gmail.com> wrote:
> Hi Lyndon,
>
> I think there is a function for that:
>
> x <- stackSelect(rb2, ind.r)
>
> I.e., for each cell, use the value in ind.r to select the layer in rb2 from
> which to take the value for the output layer.
>
> Best, Robert
>
>
>
> On Thu, Dec 1, 2011 at 5:55 PM, Lyndon Estes <lyndon.estes at gmail.com> wrote:
>>
>> 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
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>