I have two different parameters raster data, Which are temperature (X) and precipitation (Y). For temperature (X) and precipitation (Y), each parameters have 12 data files (Tiff format).

I want to perform Pearson's correlation coefficient analysis between two parameters and want output as a raster format. I have tried to do it in Arc GIS but I think in ArcGIS not providing such type computing facility.

So I am looking for a solution in an R based software package. I am new in R.

I keep my rasters in two different file, Below I have attached my folder situation which is shown below. How can i assign here path for the two different raster parameters.

enter image description here

  • 1
    I would take pause on the statistical support of this analysis. Other than addressing general curiosity, the Pearson's correlation is not well supported, statistically, for any empirical conclusions. The data are both spatially and serially autocorrelated and not necessarily multivariate normal. An alternative would be the nonparametric Kolmogorov-Smirnov test (ks.test function), which compares the equivalence of two distributions. You could make a raster of the p-values, indicating rejection of the null at a given significance. For pair-wise raster correlations look at spatialEco::hmwCorr. – Jeffrey Evans Jul 20 '17 at 16:48

I'll expose the solution with two random datasets, one for temperature and one for precipitation:


r <- raster()

t1 <- setValues(r , rnorm(n=64800,mean=16,sd=1))
t2 <- setValues(r , rnorm(n=64800,mean=18,sd=1))
t3 <- setValues(r , rnorm(n=64800,mean=20,sd=1))
t4 <- setValues(r , rnorm(n=64800,mean=22,sd=1))
t5 <- setValues(r , rnorm(n=64800,mean=24,sd=1))
t6 <- setValues(r , rnorm(n=64800,mean=26,sd=1))

p1 <- setValues(r , runif(64800,0,70))
p2 <- setValues(r , runif(64800,0,60))
p3 <- setValues(r , runif(64800,0,50))
p4 <- setValues(r , runif(64800,0,40))
p5 <- setValues(r , runif(64800,0,30))
p6 <- setValues(r , runif(64800,0,20))

temp <- stack(t1,t2,t3,t4,t5,t6)
pp <- stack(p1,p2,p3,p4,p5,p6)

Extract values for each dataset:

tempValues<- values(temp)
ppValues <- values(pp)

Create a vector to store results and make a loop for correlation computation:

corValues <- vector(mode = 'numeric')

for(i in 1:dim(tempValues)[1]){
  corValues[i] <- cor(x = tempValues[i,], y = ppValues[i,], method = 'pearson')

Create a new raster with correlation results:

correlationRaster <- setValues(r, values = corValues)



There is a unique consideration for this answer: both datasets must have the same extent, cell size and CRS.


The example was to create 6 different layers and stack them to create a multi-layer raster. In your case, just simply load each multi-layer raster with 'stack' function:

npp <- stack('D:/Test_Corr/MODIS_NPP/r001_NPP.tif')

wue <- stack('D:/Test_Corr/MODIS_WUE/r001_WUE.tif')

nppValues<- values(npp)
wueValues <- values(wue)

corValues <- vector(mode = 'numeric')

for(i in 1:dim(nppValues)[1]){
  corValues[i] <- cor(x = nppValues[i,], y = wueValues[i,], method = 'pearson')

correlationRaster <- setValues(r, values = corValues)

If datasets doesn't have the same extent, crop one of them (the bigger one) before with:

npp <- crop(npp,wue)


wue <- crop(wue,npp)
  • Thanks, i have run the code. It yields the same output which you attached. As I updated my query along with my raster files location, now how to link my raster file path with your command line t1 <- setValues(r , rnorm(n=64800,mean=16,sd=1)) ........... t12 and p1 <- setValues(r , runif(64800,0,70)) till p12. – SWAT Jul 20 '17 at 13:37
  • This seems helpful to my problem: gis.stackexchange.com/questions/277575/… Though, I am not sure how to run the correlation if I have a series of rasters. – user2543 Mar 30 '18 at 4:44
  • @user2543 is exactly what I have done here. Compare 6 rasters with other 6 rasters and compute a correlation value per pixel – aldo_tapia Mar 30 '18 at 11:43
  • Thanks. In the link I provided, I have list.files. How can I loop them from 1..200. do I stack them also? How do i do that in code? In your example, you enumated each file. It would be difficult to do that 200 times. – user2543 Mar 30 '18 at 11:46
  • @aldo_tapia: I used your code and there was an error that shows: In corValues[i] <- cor(x = redtemp[i, ], y = bluetemp[i, ... : number of items to replace is not a multiple of replacement length. My edited code can be found here: gis.stackexchange.com/questions/277575/… – user2543 Mar 31 '18 at 2:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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