R语言实现多要素偏相关分析

偏相关分析是指当两个变量同时与第三个变量相关时，将第三个变量的影响剔除，只分析另外两个变量之间相关程度的过程，判定指标是相关系数的R值。

在GIS中，偏相关分析也十分常见，我们经常需要分析某一个指数与相关环境参数的相关程度，例如NDVI与气温，降水，地形之间的相关系数。这与我们日常研究息息相关，因此掌握偏相关分析，对我们GISers比较重要。

虽然目前网络上有许多教程，但大部分是针对三个变量，且需要收费查阅，十分不便捷。

本文以标准化降水蒸散发指数SPEI与风速，降水，气温，日照时数，相对湿度等五个变量的偏相关分析为例，利用R语言进行解答。以达到大家只需更改数据存储路径，即可进行偏相关分析的目的。

话不多说，让我们开始吧！

01数据的规范

R语言进行栅格数据处理时，要求输入数据的统一规范，这包括输入栅格数据的空间参考，像元大小（分辨率），投影完全一致，我们在进行偏相关分析时，数据来源宽泛，标准不一，在数据的预处理上要特别重视，它是我们进行下一步代码操作的基础。常见的方法和工具有“投影”“重采样”“按掩膜提取”等等。这些都可以使用ArcGIS进行操作，小编在之前的文章也有过部分介绍。

02R语言代码的编写

将我们前期处理好的栅格数据按一定的命名规范，存储在本地文件夹后，我们就可以进行代码的编写了。

代码分成四部分：第一步加载所需要的包；第二步读取研究数据；第三步编写偏相关分析函数并应用；第四步导出数据。

在代码的注释中，小编已写的尽量详细，大家可以根据自己的需求来对更多或者更少变量的情况进行个性化的处理。


# 加载所需包
library(terra)
library(ppcor)

# 读取栅格数据，分别是SPEI和其他五类环境数据
MRSEI <- rast(dir("G:/demo/002/SPEI", full.names = TRUE, pattern = '.tif$'))
# 这句代码就是找到对应目录下的文件后缀名为tif……的文件，并将他们存储在MRESI这个对象中，它现在就存了9个栅格
fs <- rast(dir("G:/demo/002/五类气象数据/风速", full.names = TRUE, pattern = '.tif$'))
pre <- rast(dir("G:/demo/002/五类气象数据/降水", full.names = TRUE, pattern = '.tif$'))
tem <- rast(dir("G:/demo/002/五类气象数据/气温", full.names = TRUE, pattern = '.tif$'))
tim <- rast(dir("G:/demo/002/五类气象数据/日照时数", full.names = TRUE, pattern = '.tif$'))
shi <- rast(dir("G:/demo/002/五类气象数据/相对湿度", full.names = TRUE, pattern = '.tif$'))

# 合并处理后的栅格数据
z <- c(MRSEI, pre, tem, fs, tim, shi)#MRSEI1-9, pre：10-18,tem：19-27,fs：28-36,tim：37-45,shi：46,54

# 定义计算偏相关的函数fun_cor
fun_cor <- function(x) {
  # 如果任何一个值是 NA，则返回 NA
  if (any(is.na(x))) {
    return(c(NA, NA, NA, NA, NA))
  } else {
    # 尝试执行偏相关计算，处理可能的错误
    tryCatch({
      npre = ppcor::pcor.test(x[1:9],x[10:18],list(x[19:27],x[28:36],x[37:45],x[46:54]),method = "pearson")
      ntem = ppcor::pcor.test(x[1:9],x[19:27],list(x[10:18],x[28:36],x[37:45],x[46:54]),method = "pearson")
      nfs = ppcor::pcor.test(x[1:9],x[28:36],list(x[19:27],x[10:18],x[37:45],x[46:54]),method = "pearson")
      ntim = ppcor::pcor.test(x[1:9],x[37:45],list(x[19:27],x[28:36],x[10:18],x[46:54]),method = "pearson")
      nshi = ppcor::pcor.test(x[1:9],x[46:54],list(x[19:27],x[28:36],x[10:18],x[37:45]),method = "pearson")
      #返回降水数据的相关性和显著性检验的P值
      nprer = npre$estimate 
      nprep = npre$p.value 
      #返回气温数据的相关性和显著性检验的P值 
      ntemr = ntem$estimate
      ntemp = ntem$p.value
      nfsr = nfs$estimate
      nfsp = nfs$p.value #显著性检验的方法也有很多种
      ntimr = ntim$estimate
      ntimp = ntim$p.value
      nshir = nshi$estimate
      nship = nshi$p.value
      #在这里我就只return 五个变量的相关性了，P值就不return ,有需求的小伙伴可以自行导出
      return(c(nprer,ntemr,nfsr,ntimr,nshir))
    }, error = function(e) {
      return(c(NA, NA, NA, NA, NA))  # 如果出现错误，返回 NA
    })
  }
}
#开始运行函数
result <- terra::app(z,fun_cor,cores = 8)
#预览结果
plot(result)
将结果导出到本地，这里就导出两个为例吧
writeRaster(result[[1]],filename = "G:/demo/降水相关.tif")
writeRaster(result[[2]],filename = "G:/demo/气温相关.tif")

03 结果展示

按照本文示例代码的，我将得到五个变量的相关系数栅格图，下面是plot预览的效果。