On Wed, Jul 27, 2011 at 2:23 PM, Agustin Lobo <alobolistas@gmail.com> wrote:

> I think the problem is that I do not quite understand
> the relationship between the objects passed within overlay (s1 and s2 in
> your example) and what
> the function within overlay is actually getting.
>

It gets a matrix where the columns are layers and the rows are cells.


>
> In your example,
> x3 <- overlay(s1, s2, fun=fun3)
>
> s1 and s3 are stacks of 10x10x3 but what is fun3 getting?
> A matrix?


yes, two matrices of 100*3  (or a chunk thereof if the raster is very large)

 Note we want to calculate the median across layers,

> thus vectors of 3 elements for each pixel. As you use
> apply(x,1,median,na.rm=TRUE)
>
> it seems that x within fun3 is a matrix of 10x10 rows and 3 columns,
> is this correct?
>
> Also, why are you doing
> fun3(s1[1:3], s2[1:3])
>
>
Because I am trying it for three cells, using the type of matrix that calc
gets:  s1[1:3]


> instead of
> fun3(s1, s2)
> ?
>
>
Because fun3 may not know what to do with a RasterStack.



> Agus
>
>
Hth, Robert



> 2011/7/27 Robert J. Hijmans <r.hijmans@gmail.com>:
> > Agus,
> >
> > This problem arises because you are combing what are (for calc anyway)
> > different types of functions. The second part is called with apply, the
> > first part not. I think you can fix that using my 'fun3' (see below).
> >
> >
> > r1 <- r2 <- raster(nc=10, nr=10)
> > r1[] <- round(runif(ncell(r1))) * 248
> > r2[] <- round(runif(ncell(r2))) * 232
> >
> > s1 <- stack(r1, r1, r2)
> > s2 <- stack(r2, r2, r1)
> >
> > fun1 <- function(x,y) { x[y!=248 & y!=232] <- NA; return(x)}
> >
> > fun2 <- function(x,y) {
> > x[y!=248 & y!=232] <- NA
> > median(x,na.rm=T)
> > }
> >
> > fun3 <- function(x,y) {
> > x[y!=248 & y!=232] <- NA
> > apply(x,1,median,na.rm=TRUE)
> > }
> >
> > fun1(s1[1:3], s2[1:3])
> > fun2(s1[1:3], s2[1:3]) # not good. only returns a single value
> > fun3(s1[1:3], s2[1:3]) # looks good
> >
> > x1 <- overlay(s1, s2, fun=fun1)
> > x2 <- overlay(s1, s2, fun=fun2) # not good
> > x3 <- overlay(s1, s2, fun=fun3) # OK
> >
> >
> >
> > OK means that no error is thrown. I have not checked if the results are
> > meaningful, but I expect they are.
> >
> > Robert
> >
> >
> >
> > On Tue, Jul 26, 2011 at 8:27 AM, Agustin Lobo <alobolistas@gmail.com>
> wrote:
> >
> >> I'm trying to combine 2 operations into one to speed up a process.
> >> It's simple calculating a median raster across layers of a brick object,
> >> but
> >> the values that must be considered NA are defined by a second brick.
> >> Normally, we first set the NA in the brick and then calculate the
> median:
> >>
> >> fun1 <- function(x,y) { x[y!=248 & y!=232] <- NA; return(x)}
> >> decadaN_v2 <- overlay(decadaN, decadaS, fun=fun1, overwrite=TRUE)
> >> a <- overlay(decadaN_v2, fun=median, na.rm=T, overwrite=TRUE)
> >>
> >> But this is too slow (we actually try to do this for 36 bricks within
> >> a loop). So I tried:
> >>
> >> mimedian <- function(x,y) {
> >>  x[y!=248 & y!=232] <- NA
> >>  median(x,na.rm=T)
> >> }
> >>
> >> And then try to use overlay only once:
> >>
> >> a <- overlay(decadaN, decadaS, fun=mimedian, na.rm=T, overwrite=TRUE)
> >>
> >> but raster complains that the function is not vectorized.
> >> Is there any way around? I suspect the problem is that overlay works
> >> across layers if one brick is provided,
> >> but not if 2 bricks are provided. In other words, we want x in
> >> mimedian to be a vector of 1 pixel across layers but
> >> actually is a brick object.
> >>
> >> Agus
> >>
> >
> >        [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>

	[[alternative HTML version deleted]]

