问题引入
假设有两组栅格数据，一组代表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（），那就只能另辟蹊径了。
