I am having trouble calculating the pixel by pixel correlation coefficient between two datasets. My current code takes in two folders full of rasters and creates two independent raster stacks. These rasters all have the same cellsize and extent. I then try and take the correlation coefficient (spearman in my case) between the column values of those rasters.


r <- raster()
raster1 <- list.files(path = "data/List1", pattern = "*.tif$", full.names = T)
raster2 <- list.files(path = "data/List2", pattern = "*.tif$", full.names = T)
l1 <- stack(raster1)
l2 <- stack(raster2)

list1Values <- values(l1)
list2Values <- values(l2)
corValues <- vector(mode = 'numeric')

for (i in 1:dim(list1Values)[1]){
  corValues[i] <- cor(x = list1Values[i,], y = list2Values[i,], method = 'spearman')

corRaster <- setValues(r, values = corValues)

However at the correlation for loop it gives six error messages saying

In cor(x = list1Values[i, ], y = list2Values[i, ], method = "spearman") :
  the standard deviation is zero

Ignoring that and continuing on, the last line errors out saying

Error in setValues(r, values = corValues) : 
  length(values) is not equal to ncell(x), or to 1

My initial guess was that since the datasets have a lot of NA values (In this case it is a lot of open ocean that there is no data for), that could cause problems, so I added the

use = "complete.obs"

parameter to the cor function. This errors saying

Error in cor(x = list1Values[i, ], y = list2Values[i, ], method = "spearman",  : 
  no complete element pairs

My guess is that this last error is telling me the matrices it created don't line up, and no non-NA cells match. I have no clue how this is possible because I have been working with these rasters for years and they certainly line up. Other than that, I don't know why this isn't working.

This question is not a duplicate of Correlation/relationship between map layers in R? in any way. First I am not using point data, but instead raster. I am also interested in a simple pearson and/or spearman correlation coefficient NOT cross correlation nor spatial correlation. I'm also not using two layers, but two sets of hundreds of layers (making the statistical analysis of the correlation valid). Lastly, this code is doing fundamentally different things than that post and I am asking for specific help with the errors arising.

  • This issue was addressed just last week. In the future, please search the site before posting. Also note that you are approaching this in the same way as the previous OP, which is not correct. A correlation of two values is nonsensical. Apr 16, 2018 at 14:22
  • I will look into that issue, however note that when I did search (and I had searched for about a week prior) this particular one never showed up. Also note, this has nothing to do with a correlation of just two rasters. I have several hundred that I am correlating, and this analysis is absolutely valid.
    – maj
    Apr 16, 2018 at 18:20
  • Your code indicates that you are looking a pair-wise correlations cor(x = list1Values[i,], y = list2Values[i,]). When searching for R related content use [R] in the search term as the brackets specify the term using the letter. Apr 16, 2018 at 21:18
  • Sorry, I initially flagged the incorrect duplicate and thought that I had fixed it (apparently not). Please look at this one: gis.stackexchange.com/questions/277575/… Apr 16, 2018 at 21:20
  • Try this and just ignore the SD errors: corValues <- rep(NA,ncell(l1)); "or(i in 1:nrow(list1Values)) {corValues[i] <- cor(list1Values[i,], list2Values[i,], method = 'spearman')}"; cor.raster <- l1[[1]]; cor.raster[] <- NA; cor.raster[] <- corValues Apr 16, 2018 at 21:45

2 Answers 2


For reference, here is an approach, in memory, that uses mapply. This is one of the lesser known apply functions in R that lets you apply a function across objects.

Here, when I create objects from the raster stacks, I wrap the stack function call in list and values. This results in list objects, containing data.frames holding the raster values. This makes the data ready for mapply, which expects list objects. The use of sapply allows me to use a numeric row index 1:nrow(s1[[1]]) to aggregate the row-by-row correlations. This results in a vector of correlations, matching the rows in each lists data.frames. This vector is ordered to the cells in the source rasters and can be piped directly back into a source raster.

r <- raster(system.file("external/test.grd", package="raster"))
  s1 <- list(values( stack(r, r*0.15, r/1805.78) ))
  s2 <- list(values( stack(r^2, r/0.87*10, r/sum(values(r),na.rm=T)) ))

r[] <- as.vector(mapply(function(X,Y) { sapply(1:nrow(s1[[1]]), function(row) 
               cor(X[row,], Y[row,]))}, X=s1, Y=s2))

If I'm reading what you're trying to do correctly, a more robust approach would be to use raster::calc() on your stacked layers, as follows. This allows on-disk calculations as well as the option to speed things up with parallel processing.

# combine your two raster stacks. Now layers 1:(n/2) will be 
# correlated with layers ((n/2)+1):n. Make sure both stacks have the same 
# number of layers, i.e. nlayers(l_all) is even!
l_all <- stack(l1, l2)

# write a little function to do the correlation. The input is the vector
# of values at cell n, which has length() == nlayers(l_all)
corvec <- function(vec = NULL) {
    # 'top' half of stack
    x      = vec[1:(length(vec)/2)],
    # 'bottom' half of stack
    y      = vec[((length(vec)/2) + 1):length(vec)],
    use    = 'complete.obs',
    method = 'spearman'

corlyrs <- calc(
  fun = function(x) {
    # handle areas where all cell vals are NA, e.g. ocean
    if (all(is.na(x))) {
    } else {
      corvec(vec = x)

The above will return a rasterLayer of correlation values. It took about 30s to correlate 2 sets of 6 300x300 cell rasters on my laptop.

To speed things up with clusterR,

# keep a cpu or two spare
cpus <- max(1, (parallel::detectCores() - 1))
beginCluster(n = cpus)
corlyrs <- clusterR(x         = l_all,
                    fun       = calc,
                    args      = list(fun = function(cell) {
                      if (all(is.na(cell))) {
                      } else {
                    export    = 'corvec'

The last option may not be worth it unless your rasters are quite large. Saved me about 10 seconds on the actual calc, according to microbenchmark, but it also took about that long to establish the cluster. Benefits are much clearer with big rasters.

  • 1
    That is beautiful. I haven't messed around with the clusterR, but the base code works flawlessly. I did a pearson correlation to test the accuracy (against python code which can do this calculation except only for pearson) and the output rasters match EXACTLY. I'm fairly certain the approach Jeffery and I were eventually getting to would work as well, but this doesn't present any of the current problems. Thank you so much!
    – maj
    Apr 17, 2018 at 0:00
  • You do not have to do NA handling with the cor function. In this case NA begets NA eg., cor(c(NA,NA,NA),c(NA,NA,NA)) so, it is necessary code, albeit, good coding practice. BTW, @maj, as indicated a clear answer was on the post I linked. You can easily index the stack(s) and pipe into an empty raster. The use of multithreading is a nice addition here though. Apr 17, 2018 at 16:13
  • @orbl_soil again, this is great. I am just trying to understand exactly what it is doing. I think I understand the function itself - just getting the first half of the raster stack and correlating it with the second half, but I start to lose it when the function is actually called. Specifically after specifying the raster object to be calculated on, why is there a function(x)? Wouldn't that normally just be fun = corvec? Also, I'm trying to get a bit of feedback from the loops, by seeing what cell or something its on, but I can't find a place for it within the calculation.
    – maj
    Apr 21, 2018 at 16:48
  • @maj all good - corvec() as written will fall over on a cell where the whole stack is NA values, even though as @JeffreyEvans points out, cor() on its own looks like it should be able to handle that. More here but tbh I'm not clear on the root cause - stackoverflow.com/questions/23716967/…. I've just gotten used to adding in that little if clause to handle all-NA input vectors. I don't think there's a way to build in progress monitoring without using some kind of add-on package unfortunately.
    – obrl_soil
    Apr 22, 2018 at 22:14
  • (also the NA-handling could just as easily have been built into corvec() instead of writing the function(x) part in calc, but I forgot I needed it until the last minute)
    – obrl_soil
    Apr 22, 2018 at 22:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.