# Correlation between points and raster

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(raster)
library (gstat)
r<- raster ()

NDVI <- stack(list.files(path = userpath, pattern = 'NDVI*.*.tif',full.names = T))
proj4string(Soil)=CRS("+init=epsg:4326")
``````

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

`````` NDVItemp<-values(NDVI)
corvalues,-vector (mode='numeric')

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

correlationRaster <- setValues(NDVI[], 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 '18 at 11:23
• Hi @kamo. Where do I insert that part? – user2543 Apr 4 '18 at 12:42
• Check out the proposed solution with generated data. – Kamo Apr 4 '18 at 19:14

``````library(raster)
library(rgdal)
library(sp)

# 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. – Jeffrey Evans Apr 4 '18 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 '18 at 5:36
• Also, how do I flag significant Pearson's R? – user2543 Apr 5 '18 at 5:54
• @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 '18 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 '18 at 11:21