[R-sig-Geo] overlay with 2 bricks?

[Agustin Lobo]

ok, then, as I understand it,
given a raster stack  or brick s1 with dimensions n*c*b
(thus, 3D, being n,c,b, respectively, nb of rows, columns and layers),
overlay(s1,fun(x),...) would reshape s1 and pass it to fun as  a 2D
matrix x of dimensions
(n*c) rows and b layers,
and
given r1 as a n*c raster layer
(being n,c, respectively, nb of rows and columns),
overlay(r1,fun(x),...) would reshape r1 and pass it to fun as a vector
x of dimensions
(n*c) rows and 1 layers,

Correct?

Thus, if we want to apply an eigenvector matrix u of dimensions b*b to
s1 of dimensions (n*c*b),
we could just use

funortho <- function(x,u){ x%*%u}

a <- overlay(s1,u,fun=funortho(x,u))

as x would have dimensions (n*c) rows and b columns

Correct?

Thanks!

Agus

2011/7/27 Robert J. Hijmans <r.hijmans at gmail.com>:
> On Wed, Jul 27, 2011 at 2:23 PM, Agustin Lobo <alobolistas at 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 at 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 at 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 at r-project.org
>> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >
>>
>
>        [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>