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

[Agustin Lobo]

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.

In your example,
x3 <- overlay(s1, s2, fun=fun3)

s1 and s3 are stacks of 10x10x3 but what is fun3 getting?
A matrix? 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])

instead of
fun3(s1, s2)
?

Agus

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
>