有同学问如何使用栅格数据做冗余分析(Redundancy Analysis, RDA)，对于这种情况和类似的情况，可以参考以前分享的一篇论文：JCP|利用结构方程模型量化自然和人为因素对植被变化的影响。中国江苏省的一个案例研究，这篇论文是使用栅格和其他的多源数据做的结构方程模型，其实方法是类似的，在这里我给大家总结成了下图的标准套路：
R语言多源栅格数据模型分析标准套路
基于R语言的多源栅格数据的标准套路大约分为五个步骤：
数据输入，读取多源数据 terra包栅格裁剪、重采样、波段合成、转数据框、数据提取等：
对数据进行裁剪、重采样，由于多源数据的空间范围和分辨率可能不一致，在后续处理前最好是将空间范围和分辨率进行统一。在这里使用的是terra包，涉及函数有：裁剪 trim(mask())，重采样resample()波段合成：多源数据往往来自于多个TIF，为了方便后续的分析处理，可以使用中括号 []和c()函数进行波段选择与波段合成，将待分析的多源栅格数据组合成一个多波段的SpatRaster，方便后续的处理栅格数据转面板数据，可以使用 terra::as.data.frame()将所有的像元值转为数据框，或者extract()将像元值提取到采样点，或者zonal()分区统计等等
tidyverse来完成了，比如去除空缺值、长宽数据转换等，将数据转换为模型输入所需的可视化，对建模结果进行可视化，得到分析结果
栅格数据结构方程模型实例代码
下面给出使用栅格数据进行结构方程模型的实例代码，方法基本和前面提到的论文一致。
下面的代码符合前面的标准流程：
读取多源栅格 terra包数据提取 tidyverse数据清洗 SEM建模 可视化结构方程模型结果
需要注意的是，下面的实验结果是失败的，未达到预期目标。但是流程跑通了，把代码开放出来供大家学习参考。这种情况下就需要调整输入数据或者修改、更换模型方法，直到实验成功为止。
采样点生成可以参考下面的视频号：
#读取裁剪好的Sen+Mk结果
sggrowym_senmk_MM = rast("D:/Mongolia/NDVIClimExport/MMndviclim/Sen+MK/sggrowym_senmk_MM.tif")
tmpgrowym_senmk_MM = rast("D:/Mongolia/NDVIClimExport/MMndviclim/Sen+MK/tmpgrowym_senmk_MM.tif")
pregrowym_senmk_MM = rast("D:/Mongolia/NDVIClimExport/MMndviclim/Sen+MK/pregrowym_senmk_MM.tif")
sradgrowym_senmk_MM = rast("D:/Mongolia/NDVIClimExport/MMndviclim/Sen+MK/sradgrowym_senmk_MM.tif")
aetgrowym_senmk_MM = rast("D:/Mongolia/NDVIClimExport/MMndviclim/Sen+MK/aetgrowym_senmk_MM.tif")
soilgrowym_senmk_MM = rast("D:/Mongolia/NDVIClimExport/MMndviclim/Sen+MK/soilgrowym_senmk_MM.tif")
NTL_senmk_MM = rast("D:/Mongolia/NDVIClimExport/MMndviclim/Sen+MK/NTL_senmk_MM.tif")
#读取样本点
SP_Mo = vect("./ROI.gpkg", query = 'SELECT * FROM SP_Mo')
SP_Mo
SP_CN = vect("./ROI.gpkg", query = 'SELECT * FROM SP_CN')
SP_CN
#提取点对应的高程、起伏度、坡度、坡向值
CopDSM_Mo = terra::extract(CopDSM, SP_Mo)
CopDSM_CN = terra::extract(CopDSM, SP_CN)
#光热因子
Srad_Mo = terra::extract(sradgrowym_senmk_MM[["slope"]], SP_Mo)
Srad_CN = terra::extract(sradgrowym_senmk_MM[["slope"]], SP_CN)
Tmp_Mo = terra::extract(tmpgrowym_senmk_MM[["slope"]], SP_Mo)
Tmp_CN = terra::extract(tmpgrowym_senmk_MM[["slope"]], SP_CN)
#水分因子
Pre_Mo = terra::extract(pregrowym_senmk_MM[["slope"]], SP_Mo)
Pre_CN = terra::extract(pregrowym_senmk_MM[["slope"]], SP_CN)
Aet_Mo = terra::extract(aetgrowym_senmk_MM[["slope"]], SP_Mo)
Aet_CN = terra::extract(aetgrowym_senmk_MM[["slope"]], SP_CN)
Soil_Mo = terra::extract(soilgrowym_senmk_MM[["slope"]], SP_Mo)
Soil_CN = terra::extract(soilgrowym_senmk_MM[["slope"]], SP_CN)
#夜间灯光
NTL_Mo = terra::extract(NTL_senmk_MM[["slope"]], SP_Mo)
NTL_CN = terra::extract(NTL_senmk_MM[["slope"]], SP_CN)
#NDVI
NDVI_Mo = terra::extract(sggrowym_senmk_MM[["slope"]], SP_Mo)
NDVI_CN = terra::extract(sggrowym_senmk_MM[["slope"]], SP_CN)
#将各种参数汇总
library(tidyverse)
SEMdata_Mo = Popslope_Mo %>%
left_join(CopDSM_Mo, by = c("X"="ID")) %>%
left_join(Srad_Mo, by = c("X"="ID")) %>%
left_join(Tmp_Mo, by = c("X"="ID")) %>%
left_join(Pre_Mo, by = c("X"="ID")) %>%
left_join(Aet_Mo, by = c("X"="ID")) %>%
left_join(Soil_Mo, by = c("X"="ID")) %>%
left_join(NTL_Mo, by = c("X"="ID")) %>%
left_join(NDVI_Mo, by = c("X"="ID")) %>%
set_names(c("ID", "Lat", "Lon", "Pop", "Ele", "Relief", "Aspect", "Slope", "Srad", "Tmp", "Pre", "Aet", "Soil", "NTL", "NDVI"))
SEMdata_CN = Popslope_CN %>%
left_join(CopDSM_CN, by = c("X"="ID")) %>%
left_join(Srad_CN, by = c("X"="ID")) %>%
left_join(Tmp_CN, by = c("X"="ID")) %>%
left_join(Pre_CN, by = c("X"="ID")) %>%
left_join(Aet_CN, by = c("X"="ID")) %>%
left_join(Soil_CN, by = c("X"="ID")) %>%
left_join(NTL_CN, by = c("X"="ID")) %>%
left_join(NDVI_CN, by = c("X"="ID")) %>%
set_names(c("ID", "Lat", "Lon", "Pop", "Ele", "Relief", "Aspect", "Slope", "Srad", "Tmp", "Pre", "Aet", "Soil", "NTL", "NDVI"))
#输出结构方程模型所需数据
# write.csv(SEMdata_Mo, file = "D:/Mongolia/NDVIClimExport/SEMdata/SEMdata_Mo.csv")
# write.csv(SEMdata_CN, file = "D:/Mongolia/NDVIClimExport/SEMdata/SEMdata_CN.csv")
library(lavaan)
# 创建SEM模型语法
model <- '
# 指定潜在变量
# 定义植被变化
Vegetation =~ NDVI
#定义地理、气候、人类活动影响因子
ClimFactor =~ Pre + Aet + Soil + Srad + Tmp
GeoFactor =~ Lat + Lon + Ele + Relief + Aspect + Slope
HumanFactor =~ Pop + NTL
# 指定直接效应
Vegetation ~ ClimFactor + GeoFactor + HumanFactor
# 指定误差项
Pre ~~ Tmp
Pre ~~ Aet
Pre ~~ Soil
Pre ~~ Srad
Pop ~~ NTL
'
# 拟合模型
fit <- sem(model, data = SEMdata_CN)
summary(fit)
# 使用 fitMeasures 获取标准误差
se <- fitMeasures(fit, c("satorra.bentler.scaled", "rmsea", "srmr"))
print(se)
semPlot::semPaths(fit, "std",residuals=F,nCharNodes=4,layout="tree2",edge.label.cex = 1,sizeMan =5)
相关阅读
JCP|利用结构方程模型量化自然和人为因素对植被变化的影响。中国江苏省的一个案例研究 QGIS生成随机样点并提取栅格像元值 基于栅格数据计算重心转移模型——以NDVI重心转移为例 R语言terra包栅格数据处理基础操作与技巧详解 R语言栅格时间序列回归分析——整体和逐像元计算，并行计算，快到飞起！
更多R语言栅格数据处理课程请点击阅读原文或扫码加入地理数据科学培训班
