栅格数据实现最优参数地理探测器（OPGD）详细教程！

文章发布在下面，若对您的研究有帮助，烦请引用
K. Luo, A. Samat, T. Van de Voorde, W. Li, W. Xu and J. Abuduwaili, “Assessing Ecological Quality Dynamics and Driving Factors in the Irtysh River Basin Using AWBEI and OPGD Approaches,” in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 18, pp. 1153-1173, 2025, doi: 10.1109/JSTARS.2024.3498899.

数据准备

要探寻一堆因素对因变量的影响，首先你要确定要用哪些自变量来影响哪个因变量
想好了之后
你需要到相应的网站去下载你的研究区的自变量和因变量数据的栅格数据（可以是离散的，也可以是连续的）
后续操作是到Arcgis里对你的数据处理一下
由于不是教程的重点，这里就不展开介绍了
大致的处理就是按掩膜提取出研究区范围的栅格影像，重采样至同一分辨率
将处理后的Tif影像都存储到一个文件夹里

栅格数据读取与预处理

虽然你已经进行重采样了
但是当你用R来读取栅格影像时还有一定可能会出现错误
下面我们先来读取数据把
如果你没有这个包的话，装一下

install.packages("raster")

接下来，输入如下的代码，从而将刚刚配置好的raster包导入。

library(sp)
library(raster)

此时，在RStudio右下方的“Packages”中，可以看到raster包以及其所依赖的sp包都处于选中的状态，表明二者都已经配置成功，且完成导入
 
我们需要将存放有大量栅格图像的文件夹明确，并将其带入list.files()函数中；这一函数可以对指定路径下的文件加以遍历。其中，pattern是对文件名称加以匹配，我们用".tif$"表示只筛选出文件名称是以.tif结尾的文件；full.names表示是否将文件的全名（即路径名称加文件名称）返回，ignore.case表示是否不考虑匹配文件名称时的大小写差异。

tif_file_path <- list.files(r"(D:/lky/中科院/第三篇/驱动因素/TIFworkplace)", pattern = ".tif$", full.names = TRUE, ignore.case = TRUE)
tif_file_path
tif_file_all <- stack(tif_file_path)
tif_file_all
plot(tif_file_all)

如果没有错误的话，会这样


但如果有错误：

说明每一个tif影像的行列数不一致，这就需要我们根据所有tif影像中行列数最少的一张影像去重采样其他影像

# 定义参考文件路径
ref_file <- "D:/lky/中科院/第三篇/驱动因素/影响因子/500m/dem.tif" #这是参考影像！！！！！

# 获取所有 .tif 文件的完整路径
tif_files <- list.files("D:/lky/中科院/第三篇/驱动因素/影响因子/500m", pattern = "\\.tif$", full.names = TRUE, ignore.case = TRUE)

# 读取参考文件
ref_raster <- raster(ref_file)
ref_extent <- extent(ref_raster)

# 创建输出目录
output_dir <- "D:/lky/中科院/第三篇/驱动因素/TIFworkplace"
dir.create(output_dir, showWarnings = FALSE)

# 定义函数来调整范围并保存文件
adjust_extent <- function(file_path, ref_extent, output_dir) {
  ras <- raster(file_path)
  ras_adjusted <- projectRaster(ras, ref_raster, method = "ngb")  # 使用最近邻法调整
  
  # 构造新的文件路径，保存到指定目录下
  file_name <- tools::file_path_sans_ext(basename(file_path))
  file_path_new <- file.path(output_dir, paste0(file_name, "_adjusted.tif"))
  
  # 写入新文件
  writeRaster(ras_adjusted, file_path_new, overwrite = TRUE)
}

# 对所有 .tif 文件调用函数进行调整，并保存到指定目录
lapply(tif_files, adjust_extent, ref_extent, output_dir)

至此你的影像就没有问题了

然后再来读取一遍吧：

tif_file_path <- list.files(r"(D:/lky/中科院/第三篇/驱动因素/TIFworkplace)", pattern = ".tif$", full.names = TRUE, ignore.case = TRUE)
tif_file_path
tif_file_all <- stack(tif_file_path)
tif_file_all
plot(tif_file_all)
tif_file_all_matrix <- getValues(tif_file_all)
tif_matrix = na.omit(tif_file_all_matrix)
View(tif_matrix)
head(tif_matrix)[1:3,]

这时你的影像就已经读入R中，且转换为表格了

PS：说句题外话，你可以不看

tif_file_all_matrix <- getValues(tif_file_all)
View(tif_file_all_matrix)

这样读取后的数据是有空值的

所以tif_matrix = na.omit(tif_file_all_matrix)这句话是用来消除空值的

Optimal parameters-based geographical detector model

单因素检测

有了数据，下面就是使用GD来进行地理探测器的运行了
由于我们的数据中有连续值，但是倘若你要用geodetector包，你就需要保证输入的自变量数据已经全部为类别数据，这里也不赘述如何处理了
就是在arcgis里重分类即可
如果感兴趣可以试试用geodetector包

factor_detector("A_LCCS0", "F_LCS0", tif_frame)
factor_detector("A_LCCS0", c("DEM_Reclass", "F_LCS0"), tif_frame)

其中，"A_LCCS0"是本文中的因变量，"DEM_Reclass"与"F_LCS0"则是本文中的自变量；tif_frame则是Data Frames格式变量的名称。以上述第一句代码为例来运行，运行后稍等片刻（具体时长与数据量的大小有关），将会得到如下所示的分异及因子探测结果。
关于geodetector包的其他教程大家自己搜索吧

下面主要介绍如何使用GD包来进行
首先给大家看下我的数据

其中我的因变量是：AWBEIchange_adjusted
如果你跟我一样，数据超级多，等不及一堆代码运行，就是想一个一个的运行，那么
对于单个的类别变量你可以：

## 单个类别变量
data <- tif_matrix[, c("AWBEIchange_adjusted", "asp_adjusted")]
g1 <-gd(AWBEIchange_adjusted~asp_adjusted, data = data)  #倘若asp_adjusted是类别变量
g1

对于多个类别的加个遍历即可：

data <- tif_matrix[, c("AWBEIchange_adjusted", "asp_adjusted", "2")]
g1 <-gd(AWBEIchange_adjusted~., data = data[,1:3])  #倘若asp_adjusted是类别变量
g1

如果是单个连续遍历，则：

data <- tif_matrix[, c("AWBEIchange_adjusted", "asp_adjusted")]
discmethod <- c("equal", "natural", "quantile", "geometric", "sd")
discitv <- 5:12  # 选择5到12个区间
odc1 <- optidisc(AWBEIchange_adjusted ~ asp_adjusted, data = data, discmethod, discitv)
data$asp_adjusted <- cut(data$asp_adjusted, breaks = unique(odc1[[1]]$itv), include.lowest = TRUE)
g <- gd(AWBEIchange_adjusted ~ asp_adjusted, data = data)
g

多个的同样写个遍历就行了，这里不展示了
具体什么原理本文就不赘述了，很多文献中都有介绍，提一点的是discitv <- 5:12这里选择是的5到12中某一个分割尺度来作为最优离散程度，每一个自变量的离散程度应该是不同的，当然，如果你想可视化展示出来，可以：

data <- tif_matrix[, c("AWBEIchange_adjusted", "asp_adjusted")]
discmethod <- c("equal", "natural", "quantile", "geometric", "sd")
discitv <- 5:12  # 选择5到12个区间
odc1 <- optidisc(AWBEIchange_adjusted ~ asp_adjusted, data = data, discmethod, discitv)
dc1
plot(odc1)

可以参考下面文献作图：

届时，你已经把所有因变量与自变量的关系都运行出来了，每一对都有一个Q值和P值，根据Q值的大小可以进行后续分析，这里也不对原理做过多解释了，P值主要是用于显著性检验的。

闲话

因为作者我本人只用到单因素检测和交互检测这两种，所以文中也就只演示这两种的做法，其他的相如风险探测器: riskmean 和 gdrisk- 、 生态探测器: gdeco、空间尺度影响分析等就不展开介绍了，下面是某教程的例子，大家可以看看：





当然当然，如果你确定你的数据准确无误，时间上可以等待的话，你不妨用一步到位的方法来进行地理探测分析：

setwd("G:/OPGD")##CSV在的那个文件夹
install.packages("GD")##下载一下GD包
library(GD)##载入我们所用的GD包
data<-read.csv("FVC.csv",as.is = TRUE)##读取我们第一步做的包括了自变量和因变量的csv表
head(data)[1:3,]##就是看一眼前三行数据
discmethod <-c("equal","natural","quantile","geometric","sd")##离散化方式选择为这5中方法
discitv <-c(5:10)##分成5~10类，毕竟我们要看那种q值最大
datagdm <- gdm(FVC ~ LULC + PRE + TEM + DEM + SLOPE + POP + LIGHT + GDP,##第一个放为我们的因变量Y，其他的就是我们的自变量
              continuous_variable =c("PRE", "TEM","DEM","SLOPE","POP","LIGHT","GDP"),##这里把我们的连续变量输入在这里，类型变量就不输入
              data = data,
              discmethod = discmethod, discitv= discitv) 
datagdm
plot(datagdm)

注意了，上面代码的参数还是需要你自己改的，比如离散分割方法，尺度范围等等。当运行完代码后，所有的分析结果和图都会打印出来了。