NDVI多元线性回归残差分析，驱动因素划分，人类活动和气候变化贡献率计算及可视化（一）

内容比较多，本推文先讲数据预处理与多元线性回归残差分析计算部分，驱动因素划分，贡献率计算及可视化请关注后面的推文。点赞、分享越多，更新越快哦。点击 阅读原文 可以查看视频课程讲解。

[1] 金凯, 王飞, 韩剑桥, 等. 1982—2015年中国气候变化和人类活动对植被ndvi变化的影响[J]. 地理学报, 2020, 75(05): 961–974.

数据来源与预处理

数据来源

本复现研究使用的数据主要是GIMMS NDVI数据和气温、降水栅格数据，数据来源如下：

• GIMMS3g NDVI获取方法：

• R语言GIMMS 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数据加载、选择、查询、可视化教程合集

数据预处理

数据预处理过程总体分为三部分：

1. GIMMS NDVI数据合成与处理
2. 气温降水数据合成处理
3. 将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中合并为一个spatRaster
• tapp() 合成气温年度均值、降水年度累加值

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 
# 读取并处理每个 NetCDF 文件
for (nc_file in tmpdir) {
  # 读取 NetCDF 文件
  nc_data   
  # 裁剪到 CNroi 范围
  nc_data   nc_data   # 将裁剪后的栅格添加到列表中
  raster_list[[length(raster_list) + 1]] }

# 合并多波段的 spatRaster




    

merged_raster 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 
# 读取并处理每个 NetCDF 文件
for (nc_file in predir) {
  # 读取 NetCDF 文件
  nc_data   
  # 裁剪到 CNroi 范围
  nc_data   nc_data   # 将裁剪后的栅格添加到列表中
  raster_list[[length(raster_list) + 1]] }

# 合并多波段的 spatRaster
merged_raster 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变化的影响及相对贡献的常用方法。主要包括以下步骤：

1. 以NDVI为因变量，气温、降水为自变量，建立线性回归模型，计算模型参数
2. 以气温、降水数据和回归模型参数计算得到NDVI预测值 ，用来表示气候变化对NDVI的影响
3. 计算NDVI观测值与NDVI预测值之间的差，也就是NDVI残差 用来表示人类活动对NDVI的影响

公式如下：

式中 、 和