I edited Paulo's recommended approach to deal with NAs in the computation and it seems to work fast on a bunch of tests, including the dataset above:
stack.correlation <- function(stack1, stack2, cor.method){
# output template
cor.map <- raster(stack1)
# combine stacks
T12 <- stack(stack1,stack2)
rnlayers=nlayers(T12)
# 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, use="complete.obs"))
}
# apply the function above to each cell and write the correlation
# coefficient to the output template.
cor.map <- stackApply(T12, indices = rep(1, rnlayers),
fun = stack.sequence.cor, na.rm = FALSE)
return(cor.map)
}
cor_r=stack.correlation(stacka, stackb, "pearson")