Last active
July 16, 2023 04:57
-
-
Save hakimabdi/7308bbd6d9d94cf0a1b8 to your computer and use it in GitHub Desktop.
This function correlates time series data to produce correlation and significance grids based on the chosen method.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ############################################################################################# | |
| # 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 | |
| } | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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")