R语言：栅格序列逐像元相关性分析

原理内容：

library(terra)
library(Hmisc)
v1  = rast(array(runif(10*10*20),c(10,10,20)))#生成10*10测验栅格
v2  = rast(array(runif(10*10*20),c(10,10,20)))
z = c(v1,v2)
names(z)
nlyr(z)
fun_cor =  function(x) {
  Rs = Hmisc::rcorr(x[1:20], x[21:40], type = "spearman")
  Rx = Rs$r[2]
  Px = Rs$P[2]
  return(c(Rx, Px))
}
r_ndvi_pre = app(z, fun_cor, cores=4)
names(r_ndvi_pre) <- c('coefficient','p_value')
plot(r_ndvi_pre)

代入栅格数据并输出结果图像：

library(terra)
library(Hmisc)
ndvi <- rast(dir("F:/test/NDVI",full.names = T,pattern = '.tif$'))
tem <- rast(dir("F:/test/tem_April",full.names = T,pattern = '.tif$'))
pre <- rast(dir("F:/test/pre_April",full.names = T,pattern = '.tif$'))
z = c(ndvi,tem)
names(z)
nlyr(z)
fun_cor =  function(x) {
  Rs = Hmisc::rcorr(x[1:12], x[13:24], type = "spearman")
  Rx = Rs$r[2]
  Px = Rs$P[2]
  return(c(Rx, Px))
}
r_ndvi_pre = app(z, fun_cor, cores=4)
names(r_ndvi_pre) <- c('coefficient','p_value')
plot(r_ndvi_pre)
r <- subset(r_ndvi_pre,1)
p <- subset(r_ndvi_pre,1)
writeRaster(r, filename = "E:/test/result/r.tif")#r,输入结果路径及文件名
writeRaster(P, filename = "E:/test/result/p.tif")#p,输出结果为P值