I have a set of rasters (NDVI001.tif...NDVI200.tif) and I want to correlate these with overlain time series observed soil moisture.

This link suggests to extract point from the raster and do it graphically in Excel: I believe it can be done it R; I just do not know how.

library (gstat)
r<- raster ()

NDVI <- stack(list.files(path = userpath, pattern = 'NDVI*.*.tif',full.names = T))
Soil<-readOGR('C:\\somepath\\file.shp', 'file')

I am stuck here. I do not know how to continue. The shp has field values similar to the dates of NDVI001...200.tif.

corvalues,-vector (mode='numeric')

for (i in 1:dim(NDVI)[1]{
corValues[i] <- cor(x = NDVItemp[i,], y = Soil[i,], method = **'pearson'**)

correlationRaster <- setValues(NDVI[[1]], values = corValues)
plot(correlationRaster, xlab='latitude', ylab='longitude', main="NDVI vs Soil")  
  • Hi! I think you have to use raster::extract(NDVI, Soil) to get the NDVI values of your soil point samples. That way you will have a dataset for which correlation can be calculated column-wise.
    – Kamo
    Apr 4, 2018 at 11:23
  • Hi @kamo. Where do I insert that part?
    – user2543
    Apr 4, 2018 at 12:42
  • Check out the proposed solution with generated data.
    – Kamo
    Apr 4, 2018 at 19:14

1 Answer 1


Try this out. Follow the comments in code:


# Generate some test data
r <- raster(nrow=10,ncol=10)
values(r) <- rnorm(100)

# A raster stack with nlayers equal to number of time-steps (ti)
# In this toy dataset ti, i={1,...,5}

NDVI = stack(r,r,r,r,r)

# Now let's imagine a similar point dataset (a SpatialPointsDataframe object) with 20 points 
# also with ti, i={1,...,5}; each columns in this dataset is a different time-step

xs <- runif(20, xmin(r), xmax(r))
ys <- runif(20, ymin(r), ymax(r))

# Get some toy point data
Soil <- SpatialPointsDataFrame(coords = data.frame(xs, ys),
                       data = as.data.frame(matrix(rnorm(100), nrow=20, ncol=5)))

# Extract NDVI values in soil point data
NDVI_in_soil_samples <- extract(NDVI, Soil)

# Now let's calculate some correlations for each time-step
# Note that this solution expects that you have the same number of 
# layers in your raster stack than columns in your soil point samples

nTi = nlayers(NDVI) # Number of time-steps 

corValue <- vector(mode="numeric", length=nTi) # Correlation value here
pVal <- vector(mode="numeric", length=nTi) # p-values here

for(i in 1:nTi){

  corTest <- cor.test(Soil@data[,i], NDVI_in_soil_samples[,i], method = "pearson") 
  corValue[i] <- corTest$estimate
  pVal[i] <- corTest$p.value

To get a plot of correlation as a function of time with significant p-values (flagged as red filled points) you can do this:

pType <- c(1,16) # Point types (not-filled=1, filled=16 significant)
indPtype <- as.integer(pVal <= 0.05)+1 # set alpha of the test here (in this case alpha=0.05)
cols<- c("black","red") # Colors for points (if significant use red)

plot(1:length(corValue),corValue, type="n", xlab="Time", ylab="Pearson Correlation")
abline(h=0, lty=2, col="light grey")
lines(1:length(corValue), corValue)
points(1:length(corValue), corValue, pch=pType[indPtype], col=cols[indPtype], cex=1.5)


-- [Edited in 08/04/2018] --

If you are looking to calculate this by point (which will aggregate the time dimension) you can do the following:

nPts = nrow(Soil@data) # Number of points in the soil data

corValueByPoint <- vector(mode="numeric", length=nPts ) # Correlation value here
pValByPoint <- vector(mode="numeric", length=nPts ) # p-values here

for(i in 1:nPts ){
  # Use row data (info for one point across time)
    corTest <- cor.test(as.numeric(Soil@data[i,]), as.numeric(NDVI_in_soil_samples[i,]), method = "pearson") 

  corValueByPoint[i] <- corTest$estimate
  pValByPoint[i] <- corTest$p.value

Next, export the new soil point dataset with correlation values:

# Modify the point data to include two new columns:
# the correlation and the p-value
Soil@data <- cbind(Soil@data, 
                   data.frame(cor=corValueByPoint, pval = pValByPoint))

# Export to a shapefile
writeOGR(obj=Soil, dsn="tempdir", layer="Soil_corr", driver="ESRI Shapefile")
  • Why the loop? You could just use cor(Soil@data, NDVI_in_soil_samples) to create a correlation matrix and pairs(Soil@data, NDVI_in_soil_samples) for a scatter-plot with the correlations. Apr 4, 2018 at 21:02
  • So I was able to test the code and it worked. Thanks for that. However, I would like to get a shapefile out of the code that shows Pearson's R from the relationship of NDVI and Soil Sample. The scatterplot would be helpful, yes, but I have 52 locations of soil moisture. So plotting them in a matrix scatterplot would be confusing. But showing the Pearson's R in graduated colored symbol would be better representation.
    – user2543
    Apr 5, 2018 at 5:36
  • Also, how do I flag significant Pearson's R?
    – user2543
    Apr 5, 2018 at 5:54
  • 1
    @JeffreyEvans loop was necessary to run the cor.test() function and not cor(), that way we can get also p-values for each pairwise test.
    – Kamo
    Apr 5, 2018 at 11:12
  • Hi! @user2543 correlation works by aggregating all the points in your sample for a given time-step, so you 'lose' spatial information in this process. So, why use a shapefile if correlation represents an aggregated value for all your samples? Probably better to plot the correlation as a function of time. You can use the pVal vector to flag significant correlations. Also, if you liked the proposed solution please flag it as accepted ;-)
    – Kamo
    Apr 5, 2018 at 11:21

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.