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

[lberner at whrc.org]

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.