Correlation analysis between large number of environmental rasters

Question

I am trying to run a correlation test in R for 761 environmental rasters, but that seems to be too many to run at one time. For context, I have average monthly SST and CHL rasters over roughly 50 years that are parameters for MaxEnt species distribution modeling, and I need to know if any correlations exist between any of the raster datasets. Is there a way to load that large number of rasters for correlation analysis that I have not tried (see code below), and if not, what other options exist for determining correlation between SST and CHL rasters?

The code I've tried:

library(ENMTools)

library(raster)

setwd(“environmental rasters location”)

rastlist <- list.files(path = “file location”, pattern = ‘.asc’, all.files = TRUE, full.names = FALSE)

allrasters <- lapply(rastlist, raster)

allrasters_stacked <- stack(allrasters)

correlations <- raster.cor.matrix(allrasters_stacked)

This didn't work so I tried:

r1 <- raster("Chl_Apr2003.1km.asc")

correlations <- cor(sampleRandom(allrasters_stacked, size = ncell(r1)*0.05), method = "spearman")

and got the error:

Cannot allocate vector of size 7.6 Gb

Answer 1

If you have enough memory to read in and correlate two rasters you can do it one pair at a time.

I've created a sample set of four .asc raster files.

> rastlist
[1] "r1.asc" "r2.asc" "r3.asc" "r4.asc"

Processing as per your script I can do them all in one go to get:

> raster.cor.matrix(allrasters_stacked)
             r1           r2           r3          r4
r1  1.000000000  0.006363855  0.002696825 -0.02095291
r2  0.006363855  1.000000000 -0.031109053 -0.25001460
r3  0.002696825 -0.031109053  1.000000000 -0.17824323
r4 -0.020952905 -0.250014600 -0.178243229  1.00000000

But I can also loop over the off-diagonal terms, creating two rasters from the file name, and computing each element separately:

res = lapply(1:(length(rastlist)-1), function(i){
    ri = raster(rastlist[i])
    lapply((i+1):length(rastlist), function(j){
        c(i, j, 
        raster.cor.matrix(stack(ri, raster(rastlist[j])))[1,2])
    })
})

matrix(unlist(res), ncol=3, byrow=TRUE)

Giving:

     [,1] [,2]         [,3]
[1,]    1    2  0.006363855
[2,]    1    3  0.002696825
[3,]    1    4 -0.020952905
[4,]    2    3 -0.031109053
[5,]    2    4 -0.250014600
[6,]    3    4 -0.178243229

which are the same numbers as the correlation matrix only in row, column, value form.

Your problem might then be time - for 761 rasters you will be computing 289180 (761*760/2) correlation pairs. But the process is trivially parallelizable so if you have a cluster handy you can spread the load.

Depending on your memory available you could possibly do the correlation matrix in blocks, or by rows or columns with a similar looping approach.

Sidenotes: You should look at using the terra package instead of raster which might be faster and lighter, and also using a better raster file format (such as GeoTIFF) than .asc files, but if that's what you're given then that's what you have to work with, although my looping approach will read in each file every time meaning a more efficient raster file format might give massive speedups - run through your .asc files and save as .tiff and see what happens. Try small examples (like I've done) to assess correctness and see how the problem scales before wasting CPU time on buggy code or code that might take too long...