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 <[email protected]> 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 > [email protected] > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
