Skip to content

overlay with 2 bricks?

6 messages · Agustin Lobo, Robert J. Hijmans

Original
Agustin Lobo
# 
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
1 day later
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>:
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>: