Skip to content

questions on RasterStack/Brick

3 messages · lberner at whrc.org, Robert J. Hijmans

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.
3 days later
lberner at whrc.org
# 
Apologies-- there was an error in the stack correlation function of my last
post. Below is an updated version of the function.  

stack.correlation <- function(stack1, stack2, cor.method, na.val, na.rm =
T){
  # output template
  cor.map <- raster(stack1)
  # combine stacks
  combined.stack <- stack(stack1,stack2)
  NAvalue(combined.stack) = na.val
  
  # The function takes a vector of cell values (z-dimension), partitions it
in half, then correlates
  # the two sections, returning the correlation coefficient.
  stack.sequence.cor <- function(full.vec,na.rm=T){
    vec1<-full.vec[1:(length(full.vec)/2)]
    vec2<-full.vec[(length(full.vec)/2+1):length(full.vec)]
    return(cor(vec1,vec2, method=cor.method))
    }

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

example call: 
my.cor.map <- stack.correlation(august_temp, july_temp, "spearman", -9999)
#plot(my.cor.map)


Logan


--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/questions-on-RasterStack-Brick-tp5553580p6634917.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
Robert J. Hijmans
# 
Thanks Logan, 
That is very useful, but I would not use stackApply, but the simpler calc.
Below a reworked version of your function:


stackcor <- function(s1, s2, method='spearman') {
	mycor <- function(v) {
		x <- v[1:split]
		y <- v[(split+1):(2*split)]
		cor(x, y, method=method)
	}
	s <- stack(s1, s2)
	split <- nlayers(s)/2
	calc(s, fun=mycor )
}

Best, Robert

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