本文的目标是复现下面《地理学报》的残差分析部分，主要是两幅图，驱动因素的空间分布图，和人类活动、气候变化贡献率的分布。由于使用的数据略有差异，结果可能有所不同。
内容比较多，本推文先讲数据预处理与多元线性回归残差分析计算部分，驱动因素划分，贡献率计算及可视化请关注后面的推文。点赞、分享越多，更新越快哦。点击阅读原文可以查看视频课程讲解。
[1] 金凯, 王飞, 韩剑桥, 等. 1982—2015年中国气候变化和人类活动对植被ndvi变化的影响[J]. 地理学报, 2020, 75(05): 961–974.
数据来源与预处理
数据来源
本复现研究使用的数据主要是GIMMS NDVI数据和气温、降水栅格数据，数据来源如下：
GIMMS3g NDVI获取方法：
1901-2022年中国1km分辨率逐月平均气温数据集
http://www.geodata.cn/data/datadetails.html?dataguid=164304785536614&docId=686 1901-2022年中国1km分辨率逐月降水量数据集
http://www.geodata.cn/data/datadetails.html?dataguid=192891852410344&docId=687 气温、降水数据获取请参阅：如何找到靠谱的科研数据？免费的国家数据中心为你服务！
中国陆地面要素矢量：使用国家1：100万基础地理数据库中的数据合成得到，在QGIS里面另存为一个gpkg，方便调用，具体参阅：
2021版国家基础地理数据库批量下载与拼接方法 1：100万基础地理数据库怎么看？怎么用？能干嘛？ QGIS数据加载、选择、查询、可视化教程合集
数据预处理
数据预处理过程总体分为三部分：
GIMMS NDVI数据合成与处理 气温降水数据合成处理 将GIMMS NDVI和气温降水数据时序、分辨率统一起来，作为后续多元线性回归和残差分析的输入
GIMMS NDVI数据合成与处理
GIMMS NDVI的处理在前面写过，不再多解释，下面直接给出代码，不明白的可以扫码看课程讲解。
library(gimms)
# 下载GIMMS数据
downloadGimms(x=1981,y=2015,dsn="./GIMMSDOWNLOAD/")
#获取GIMMS数据列表
gimms_filenames = list.files(path = "./GIMMSDOWNLOAD", pattern = '.nc4')
gimms_files = paste0("./GIMMSDOWNLOAD/", gimms_filenames)
#读取研究区范围，使用WGS84坐标系的SHP
library(sf)
cmrshp = read_sf(dsn="D:/TEMP/TEMPROI.gpkg", query = 'SELECT * FROM CNpolygon')
cmrshp = as_Spatial(cmrshp)
#读取GIMMS数据，栅格化
gimms_rasterq = rasterizeGimms(x=gimms_files, ext = cmrshp, keep = 0, cores = 4)
mvc = monthlyComposite(gimms_rasterq, indices=monthlyIndices(gimms_files), cores=28)
yearIndices = sort(rep(1:34, 12))
ymc = monthlyComposite(mvc, fun=mean, indices=yearIndices, cores=28)
ymcnames = seq(as.Date("1982-01-01"), as.Date("2015-12-31"), "years")
ymc
plot(ymc)
library(terra)
ymc = rast(ymc)
names(ymc) = ymcnames
ymccn = trim(mask(ymc, CNroi))
plot(ymccn)
writeRaster(ymccn, filename = paste0("./GIMMSYM/", ymccn@ptr[["names"]], ".tif"))
气温降水数据处理
下面写一下气温、降水数据的处理。将从国家地球系统科学数据中心下载的气温、降水数据分别放到一个文件夹中，解压缩后如下图所示，就是下面代码的输入。
下面的代码主要解决了以下几个问题，气温、降水的处理方法是一致的，同样类型的代码写了两遍，变量名不一致而已。 fs::dir_ls()批量读取nc文件的文件路径，用来批量读取nc数据使用 for循环批量读取nc，在这里进行了数据的裁剪、重采样，重采样是为了解决由于彭老师提供的数据范围不同年份有所差异的问题，只有统一了数据的范围才能将所有的nc读取后在R中合并为一个spatRastertapp()合成气温年度均值、降水年度累加值
library(terra)
CNroi = vect("D:/TEMP/TEMPROI.gpkg", query = 'SELECT * FROM CNpolygon')
plot(CNroi)
#读取气温nc文件夹
tmpdir = fs::dir_ls(path = "G:/Geodata/PengClim/tmp/", regexp = ".nc$")
tmpdir
tmp1 = rast(tmpdir[1]) #读取第一个nc数据
tmp1 #查看数据信息
plot(tmp1) #预览数据
# 创建一个空的多波段 spatRaster 列表
raster_list <- list()
# 读取并处理每个 NetCDF 文件
for (nc_file in tmpdir) {
# 读取 NetCDF 文件
nc_data <- rast(nc_file)
# 裁剪到 CNroi 范围
nc_data <- crop(nc_data, CNroi)
nc_data <- terra::resample(nc_data, rast(tmpdir[1]))
# 将裁剪后的栅格添加到列表中
raster_list[[length(raster_list) + 1]] <- nc_data
}
# 合并多波段的 spatRaster
merged_raster <- rast(raster_list)
tsname = seq(as.Date("1982-01-01"), as.Date("2022-12-31"), "month")
tsname
names(merged_raster) = tsname
# 打印合并后的 spatRaster 信息
merged_raster
#计算温度年度均值
yearIndices = sort(rep(1:41, 12))
yearIndices
tmpYearmean = tapp(merged_raster, index = yearIndices, fun = mean, cores = 28)
names(tmpYearmean) = seq(1982, 2022, 1)
tmpYearmean = tmpYearmean/10 #将0.1℃的单位转换为℃
tmpYearmean
plot(tmpYearmean)
writeRaster(tmpYearmean, filename = paste0("./PengClim/tmpYearmean/", tmpYearmean@ptr[["names"]], ".tif"))
#读取降水栅格
predir = fs::dir_ls(path = "G:/Geodata/PengClim/pre/", regexp = ".nc$")
# 创建一个空的多波段 spatRaster 列表
raster_list <- list()
# 读取并处理每个 NetCDF 文件
for (nc_file in predir) {
# 读取 NetCDF 文件
nc_data <- rast(nc_file)
# 裁剪到 CNroi 范围
nc_data <- crop(nc_data, CNroi)
nc_data <- terra::resample(nc_data, rast(predir[1]))
# 将裁剪后的栅格添加到列表中
raster_list[[length(raster_list) + 1]] <- nc_data
}
# 合并多波段的 spatRaster
merged_raster <- rast(raster_list)
tsname = seq(as.Date("1982-01-01"), as.Date("2022-12-31"), "month")
tsname
names(merged_raster) = tsname
# 打印合并后的 spatRaster 信息
merged_raster
#计算温度年度均值
yearIndices = sort(rep(1:41, 12))
yearIndices
preYearsum = tapp(merged_raster, index = yearIndices, fun = sum, cores = 28)
names(preYearsum) = seq(1982, 2022, 1)
preYearsum = preYearsum/10 #将0.1mm单位改为1mm
preYearsum
plot(preYearsum)
writeRaster(preYearsum, filename = paste0("./PengClim/preYearsum/", preYearsum@ptr[["names"]], ".tif"))
需要注意的是，上面的气温、降水数据我是把1982-2022年也就是写这个推文的时候所有的1982年后的气温、降水数据都合并了，数据量巨大，合成过程会占用非常高的内存，最高使用到了近256GB，如果没有这么大内存的电脑，谨慎使用，或者分批进行合成，不要一口气把所有的数据都进行处理。
同时间序列数据选取与重采样
由于GIMMS NDVI的分辨率大约为8km，气温降水约1km，而且时间序列不同，需要对数据进行一些波段的选取和重采样操作。
#读取气温、降水年度数据
tmpYearmean = rast(fs::dir_ls("./PengClim/tmpYearmean/", regexp = ".tif$")[1:34])
preYearsum = rast(fs::dir_ls("./PengClim/preYearsum/", regexp = ".tif$")[1:34])
tmpYearmean
preYearsum
ymccn
#将气温、降水重采样为和GIMMS NDVI一致的分辨率
tmpYearmean = terra::resample(tmpYearmean, ymccn)
preYearsum = terra::resample(preYearsum, ymccn)
多元线性回归与残差分析
计算原理
多元回归残差分析是研究人类活动和气候变化对NDVI变化的影响及相对贡献的常用方法。主要包括以下步骤：
以NDVI为因变量，气温、降水为自变量，建立线性回归模型，计算模型参数 以气温、降水数据和回归模型参数计算得到NDVI预测值，用来表示气候变化对NDVI的影响 计算NDVI观测值与NDVI预测值之间的差，也就是NDVI残差用来表示人类活动对NDVI的影响
公式如下：
式中、和分别代表NDVI预测值、NDVI残差值和NDVI观测值，a、b、c分别对应线性回归模型的气温、降水对应的系数和模型截距
具体实现
实现思路还是栅格计算，逐像元计算多元线性回归，后面的残差分析其实就是最基本的四则运算了。
#多元线性回归计算系数
fun_mlr = function(x){
if(length(na.omit(x))<102) return(c(NA, NA, NA, NA)) #删除数据不连续含有NA的像元，NA数量必须和后面return的变量数一致
lm1 = lm(x[1:34]~x[35:68]+x[69:102])
a1 = summary(lm1)
a3 = a1$fstatistic
pvalue = 1-pf(a3[1], a3[2], a3[3])
intercept = a1$coefficients[1]
a_tmp = a1$coefficients[2]
b_pre = a1$coefficients[3]
return(c(intercept, a_tmp, b_pre, pvalue))
}
#输入待计算的数据
#3个参数,NDVI，气温、降水
lmrast = c(ymccn, tmpYearmean, preYearsum)
plot(lmrast) #预览待线性回归计算的数据
lmrast
writeRaster(lmrast, filename = "./lmrast.tif")
rast_mlr = app(lmrast, fun_mlr, cores = 56) #多元线性回归计算
names(rast_mlr) = c("intercept", "a_tmp", "b_pre", "pvalue")
plot(rast_mlr)
summary(rast_mlr)
以上的代码，我们就实现了逐像元多元线性回归的计算，计算结果预览如下：
intercept截距，对应公式中的c a_tmp，对应公式中的气温系数a b_pre，对应公式中的降水系数b
#计算NDVI预测值
ndvi_pred = rast_mlr[["a_tmp"]]*tmpYearmean+rast_mlr[["b_pre"]]*preYearsum+rast_mlr[["intercept"]]
names(ndvi_pred) = seq(1982,2015,1)
ndvi_pred
plot(ndvi_pred)
#计算NDVI残差
ndvi_res = ymccn-ndvi_pred
ndvi_res
以上就完成了NDVI多元线性回归残差分析的计算，但是距离论文中的结果图还有一定的距离，后面还需进一步的处理。
可能出现的问题
栅格计算可能会出现如下的报错：
错误: Not compatible with requested type: [type=list; target=double]. 此外: Warning message: In matrix(r, ncol = nlyr(out), byrow = TRUE) : 数据长度[385293]不是矩阵行数[96324]的整倍
参考文献
金凯, 王飞, 韩剑桥, 等. 1982—2015年中国气候变化和人类活动对植被ndvi变化的影响[J]. 地理学报, 2020, 75(05): 961–974. https://github.com/rspatial/terra/issues/962
