问题引入

假设有两组栅格数据，一组代表2019年中国每月降雨量，一组代表2019年中国每月植被叶面积指数（LAI）。想要得到中国月降水量与LAI的相关性分布，那么需要对两组栅格数据对应的栅格点进行逐栅格的相关性分析。

R包加载

library(raster)
library(broom)
library(rgdal)

数据导入

将降水数据导入栅格栈中，这个过程可以理解为将降水数据按时间顺序从上到下堆叠。同理，按相同的时间顺序将LAI数据堆叠。值得一提的是，stack（）函数在堆叠栅格数据时是按文件名拼音和数字大小顺序自动堆叠的，具体规则可以亲自尝试。最后，将这两个栅格栈合并成一个。

setwd("D:/1data/rain")                                    #降水数据所在文件夹
raster1<-stack(list.files(pattern='*.tif$'))         #堆栈降水栅格数据
setwd("D:/1data/LAI")                                    #LAI数据所在文件夹
raster1<-stack(list.files(pattern='*.tif$'))         #堆栈LAI栅格数据
r_all<-stack(raster1,raster2)

编写函数

对相关性分析函数稍作改变。

corvec <- function(vec = NULL) {
  cor(
    # 上半部分栅格栈
    x      = vec[1:(length(vec)/2)],
    # 下半部分栅格栈
    y      = vec[((length(vec)/2) + 1):length(vec)],
    use    = 'complete.obs',
    #如果是线性相关性分析改成‘person’
    method = 'spearman'           
  )
}

逐栅格计算

corlyrs <- calc(
  r_all,
  fun = function(x) {
    # 对缺失值进行处理
    if (all(is.na(x))) {
      NA_real_
    } else {
      corvec(vec = x)
    }
  }
)

栅格输出

writeRaster(corlyrs,filename = file.path( "E:/1data","降水与LAI相关性分析.tif"),format="GTiff", overwrite=TRUE)

小结

以上方法是可以推广的，线性回归函数lm（）和相关性分析函数cor（）的输入都可以是向量，因此只要函数支持向量输入，理论上讲都可以类比上述过程实现。但是如果函数只支持数据框输入，如gbm包中的函数gbm（），那就只能另辟蹊径了。