Dear all Can I use residual cokriging for some variables?? If yes,how?? On 12/2/11, Robert J. Hijmans <[email protected]> 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 <[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 >
-- Saman Monfared MSc. Student of statistics, Department of science, Shiraz university Shiraz-Iran _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
