# [R-sig-Geo] questions on RasterStack/Brick

lberner at whrc.org lberner at whrc.org
Tue Jul 26 22:21:16 CEST 2011

```Dear all,

Robert suggested that, with a little patience, the following code could be
used to correlate two raster stacks:

r <- raster(s1)
for (i in 1:ncell(r)) {
r[i] <- cor(as.vector(cellValues(s1, i)), as.vector(cellValues(s2, i)))
}

Not having that much patience, my co-worked Pieter Beck and I wrote a
function to correlate two raster stacks of same x, y, and z extents. The
function is considerably faster than the code above:

stack.correlation <- function(stack1, stack2, cor.method, na.val, na.rm =
T){
# output template
cor.map <- raster(stack1)
# combine stacks
T12 <- stack(stack1,stack2)
NAvalue(T12) = na.val

# the function takes a vector, partitions it in half, then correlates
# the two sections, returning the correlation coefficient.
stack.sequence.cor <- function(myvec,na.rm=T){
myvecT1<-myvec[1:(length(myvec)/2)]
myvecT2<-myvec[(length(myvec)/2+1):length(myvec)]
return(cor(myvecT1,myvecT2, method =  cor.method))
}

# apply the function above to each cell and write the correlation
# coefficient to the output template.
cor.map <- stackApply(T12, indices = rep(1, length(jul_mckenney)),
fun = stack.sequence.cor, na.rm = TRUE)

return(cor.map)
}

example call:
my.cor.map <- stack.correlation(tmp.jul, tmp.aug, cor.method = "spearman",
na.val = -9999)

Hope this helps,

Logan Berner
lberner at whrc.org
The Woods Hole Research Center

--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/questions-on-RasterStack-Brick-tp5553580p6623580.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

```