I have six raster datasets with matching resolution and extent.

I want to analyze the relation between those rasters by calculating pairwise Spearman's correlation coefficients for all combinations of two raster layers.

The aim is to output a correlation matrix that shows the coefficient for each of the combinations of bands.

I tried this https://www.rdocumentation.org/packages/raster/versions/2.9-23/topics/layerStats, it worked, but not for spearman.

myfolder<- "G:/name/name_2/test"
r_path <- file.path(myfolder, grep(".tif$",
                                              all.files = F),
                                   ignore.case = TRUE, value = TRUE))
mystack <- raster::stack(r_path) #https://www.researchgate.net/post/How_do_I_make_a_raster_correlation_map_using_different_raster_response_ecosystem_service_value_and_explanatory_temperature_NDVI_raster_data
raster::layerStats(mystack, 'pearson', na.rm=T)

I found the function corLocal (Spearman correlation between two rasters in R), but it performed a spearman correlation between two raster datasets and I need a correlation matrix for six raster datasets.

temp = list.files(pattern="*.tif$") #source https://stackoverflow.com/questions/52746936/how-to-efficiently-import-multiple-raster-tif-files-into-r
corLocal(raster1, raster2, method="spearman")
Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘corLocal’ for signature ‘"character", "character"’

Are there any suggestions?

closed as off-topic by Spacedman, MrXsquared, TomazicM, PolyGeo Oct 4 at 20:03

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "When seeking help to debug/write/improve code always provide the desired behavior, a specific problem/error and the shortest code (as formatted text, not pictures) needed to reproduce it in the question body. Providing a clear problem statement and a code attempt helps others to help you." – PolyGeo
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 1
    What exactly are you trying to correlate across the six rasters? Correlations (like Spearman) need two sets of observations. – Spacedman Oct 4 at 16:30

First convert the rasters to variables in the same dataframe, then calculate the pairwise correlations and use the package 'corrplot' to display the results in a matrix.


#read in rasters
r1 <- raster("IMG_0003_1.tif")
r2 <- raster("IMG_0003_2.tif")
r3 <- raster("IMG_0003_3.tif")
r4 <- raster("IMG_0003_4.tif")
r5 <- raster("IMG_0003_5.tif")

#stack raster layers
st <- stack(r1, r2, r3, r4, r5)

#subsample 5% of pixels and calculate pairwise correlations
cor<- cor(sampleRandom(st, size= ncell(r1) * 0.05 ), method = "spearman")

#plot correlation matrix
df <- corrplot(cor, method = "number")


Edited to analyze subsample of pixels rather than entire population as recommended by @JefferyEvans.

  • 2
    There is a good reason that the layerStats function uses a sample. A raster would be considered a population and has parametric convergence characteristics. As such, a correlation based on the entire population would be biased towards a Gaussian process regardless of any inherent distributional characteristics of the data. The curse of large numbers! With this approach you can functionally say “at least we know the correlation is not 0”. You could easily add sampling to this method using an index i <- sample(1:ncell(r), 10000) then rs <- r[i] recycling the index across all the rasters. – Jeffrey Evans Oct 4 at 20:53
  • Thanks. I added the sampling to my answer. – Cory G. Oct 4 at 21:06
  • How about simply stacking the rasters and using: cor(sampleRandom(r, size= ncell(r) * 0.05 )) – Jeffrey Evans Oct 4 at 22:01
  • @JeffreyEvans Much cleaner that way and the sample size has context now. – Cory G. Oct 6 at 0:38
  • @JeffreyEvans, Cory G. thank you very much, the code works perfectly. – nora Oct 7 at 11:34

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