This function correlates time series data to produce correlation and significance grids based on the chosen method.

gridcorts

#############################################################################################
# title         : Gridded Correlation of Time Series Raster Data (Gridcorts)
# purpose       : Pixelwise time series correlation and significance based on Pearson's r, 
#               : Spearman's rho or Kendall's tau
# author        : Abdulhakim Abdi (@HakimAbdi)
# input         : Raster brick comprising two time series of equal number of layers
# output        : Raster object of pixelwise correlation, significance or both (i.e. brick) 
#               : based on the chosen method
# update        : Minor based suggestion from Tao Xiong of Northeast Normal University in China.
# data          : Test data: https://1drv.ms/u/s!AsHsKb_qtbkwgvoQuHnPaFazPr_XnA?e=fJ9fue
#############################################################################################
# Please cite this paper if you use this code in a publication
#############################################################################################
# CITATION: Abdi, Abdulhakim M., et al. "The El Niño–La Niña cycle and recent trends in supply 
# and demand of net primary productivity in African drylands." Climatic Change 138.1-2 (2016): 111-125. 
# https://link.springer.com/article/10.1007/s10584-016-1730-1
#############################################################################################

gridcorts <- function(rasterstack, method, type=c("corel","pval","both")){
  # Values for (layers, ncell, ncol, nrow, method, crs, extent) come straight from the input raster stack
  # e.g. nlayers(rasterstack), ncell(rasterstack)... etc.
  print(paste("Start Gridcorts:",Sys.time()))
  print("Loading parameters")
  layers=nlayers(rasterstack);ncell=ncell(rasterstack);
  ncol=ncol(rasterstack);nrow=nrow(rasterstack);crs=crs(rasterstack);
  extent=extent(rasterstack);pb = txtProgressBar(min = 0, max = ncell, initial = 0)
  print("Done loading parameters")
  mtrx <- as.matrix(rasterstack,ncol=layers)
  empt <- matrix(nrow=ncell, ncol=2)
  print("Initiating loop operation")
  if (type == "corel"){
    for (i in 1:ncell){
      setTxtProgressBar(pb,i)
      if (all(is.na(mtrx[i,1:(layers/2)])) | all(is.na(mtrx[i,((layers/2)+1):layers]))){ 
        empt[i,1] <- NA 
      } else 
        if (sum(!is.na(mtrx[i,1:(layers/2)]/mtrx[i,((layers/2)+1):layers])) < 4 ){
          empt[i,1] <- NA 
        } else 
          empt[i,1] <- as.numeric(cor.test(mtrx[i,1:(layers/2)], mtrx[i,((layers/2)+1):layers],method=method)$estimate)
    }
    print("Creating empty raster")
    corel <- raster(nrows=nrow,ncols=ncol,crs=crs)
    extent(corel) <- extent
    print("Populating correlation raster")
    values(corel) <- empt[,1]
    print(paste("Ending Gridcorts on",Sys.time()))
    corel
  } 
  else
    if (type == "pval"){
      for (i in 1:ncell){
        setTxtProgressBar(pb,i)
        if (all(is.na(mtrx[i,1:(layers/2)])) | all(is.na(mtrx[i,((layers/2)+1):layers]))){ 
          empt[i,2] <- NA 
        } else 
          if (sum(!is.na(mtrx[i,1:(layers/2)]/mtrx[i,((layers/2)+1):layers])) < 4 ){
            empt[i,2] <- NA 
          } else 
            empt[i,2] <- as.numeric(cor.test(mtrx[i,1:(layers/2)], mtrx[i,((layers/2)+1):layers],method=method)$p.value)
      }
      pval <- raster(nrows=nrow,ncols=ncol,crs=crs)
      extent(pval) <- extent
      print("Populating significance raster")
      values(pval) <- empt[,2]
      print(paste("Ending Gridcorts on",Sys.time()))
      pval
    }
  else
    if (type == "both"){
      for (i in 1:ncell){
        setTxtProgressBar(pb,i)
        if (all(is.na(mtrx[i,1:(layers/2)])) | all(is.na(mtrx[i,((layers/2)+1):layers]))){ 
          empt[i,] <- NA 
        } else 
          if (sum(!is.na(mtrx[i,1:(layers/2)]/mtrx[i,((layers/2)+1):layers])) < 4 ){
            empt[i,] <- NA 
          } else {
            empt[i,1] <- as.numeric(cor.test(mtrx[i,1:(layers/2)], mtrx[i,((layers/2)+1):layers],method=method)$estimate) 
            empt[i,2] <- as.numeric(cor.test(mtrx[i,1:(layers/2)], mtrx[i,((layers/2)+1):layers],method=method)$p.value)
          }
      }
      c <- raster(nrows=nrow,ncols=ncol,crs=crs)
      p <- raster(nrows=nrow,ncols=ncol,crs=crs)
      print("Populating raster brick")
      values(c) <- empt[,1]
      values(p) <- empt[,2]
      brk <- brick(c,p)
      extent(brk) <- extent
      names(brk) <- c("Correlation","Pvalue")
      print(paste("Ending Gridcorts on",Sys.time()))
      brk
    }
}

Thank you very much. Please tell me how to get the correlation image?
require(raster)
npp<-raster("C:/Users/User/OneDrive/Desktop/Corr/npp.tif")
rainfall<-raster("C:/Users/User/OneDrive/Desktop/Corr/rain.tif")
rstack <- stack(npp,rainfall)
gridcorts <- function(rasterstack, method, type=c("corel","pval","both")){

Values for (layers, ncell, ncol, nrow, method, crs, extent) come straight from the input raster stack

e.g. nlayers(rasterstack), ncell(rasterstack)... etc.

print(paste("Start Gridcorts:",Sys.time()))
print("Loading parameters")
layers=nlayers(rasterstack);ncell=ncell(rasterstack);
ncol=ncol(rasterstack);nrow=nrow(rasterstack);crs=crs(rasterstack);
extent=extent(rasterstack);pb = txtProgressBar(min = 0, max = ncell, initial = 0)
print("Done loading parameters")
mtrx <- as.matrix(rasterstack,ncol=layers)
empt <- matrix(nrow=ncell, ncol=2)
print("Initiating loop operation")
if (type == "corel"){
for (i in 1:ncell){
setTxtProgressBar(pb,i)
if (all(is.na(mtrx[i,1:(layers/2)])) | all(is.na(mtrx[i,((layers/2)+1):layers]))){
empt[i,1] <- NA
} else
if (sum(!is.na(mtrx[i,1:(layers/2)]/mtrx[i,((layers/2)+1):layers])) < 4 ){
empt[i,1] <- NA
} else
empt[i,1] <- as.numeric(cor.test(mtrx[i,1:(layers/2)], mtrx[i,((layers/2)+1):layers],method=method)$estimate)
}
print("Creating empty raster")
corel <- raster(nrows=nrow,ncols=ncol,crs=crs)
extent(corel) <- extent
print("Populating correlation raster")
values(corel) <- empt[,1]
print(paste("Ending Gridcorts on",Sys.time()))
corel
}
else
if (type == "pval"){
for (i in 1:ncell){
setTxtProgressBar(pb,i)
if (all(is.na(mtrx[i,1:(layers/2)])) | all(is.na(mtrx[i,((layers/2)+1):layers]))){
empt[i,2] <- NA
} else
if (sum(!is.na(mtrx[i,1:(layers/2)]/mtrx[i,((layers/2)+1):layers])) < 4 ){
empt[i,2] <- NA
} else
empt[i,2] <- as.numeric(cor.test(mtrx[i,1:(layers/2)], mtrx[i,((layers/2)+1):layers],method=method)$p.value)
}
pval <- raster(nrows=nrow,ncols=ncol,crs=crs)
extent(pval) <- extent
print("Populating significance raster")
values(pval) <- empt[,2]
print(paste("Ending Gridcorts on",Sys.time()))
pval
}
else
if (type == "both"){
for (i in 1:ncell){
setTxtProgressBar(pb,i)
if (all(is.na(mtrx[i,1:(layers/2)])) | all(is.na(mtrx[i,((layers/2)+1):layers]))){
empt[i,] <- NA
} else
if (sum(!is.na(mtrx[i,1:(layers/2)]/mtrx[i,((layers/2)+1):layers])) < 4 ){
empt[i,] <- NA
} else {
empt[i,1] <- as.numeric(cor.test(mtrx[i,1:(layers/2)], mtrx[i,((layers/2)+1):layers],method=method)$estimate)
empt[i,2] <- as.numeric(cor.test(mtrx[i,1:(layers/2)], mtrx[i,((layers/2)+1):layers],method=method)$p.value)
}
}
c <- raster(nrows=nrow,ncols=ncol,crs=crs)
p <- raster(nrows=nrow,ncols=ncol,crs=crs)
print("Populating raster brick")
values(c) <- empt[,1]
values(p) <- empt[,2]
brk <- brick(c,p)
extent(brk) <- extent
names(brk) <- c("Correlation","Pvalue")
print(paste("Ending Gridcorts on",Sys.time()))
brk
}
}

cg1 <- gridcorts(rasterstack = rstack, method = "spearman",type="both")