地理探测器R语言实现：geodetector

2024-03-18 10:28 由疯狂学习GIS 发表于 #后端开发

本文介绍基于R语言中的geodetector包，依据多张栅格图像数据，实现地理探测器（Geodetector）操作的详细方法。

需要说明的是，在R语言中进行地理探测器操作，可以分别通过geodetector包、GD包等2个包实现。其中，geodetector包是地理探测器模型的原作者团队早先开发的，其需要保证输入的自变量数据已经全部为类别数据；而GD包则是另外一位学者开发的，其可以自动实现自变量数据的最优离散化方法选取与执行——即我们可以直接把自变量带入这一包中，无需额外进行数据的离散化。本文介绍的是基于前者，即geodetector包实现地理探测器的具体操作；基于后者的方法，我们将在后期的博客中介绍。此外，如果希望基于Excel实现地理探测器，大家可以参考地理探测器Geodetector下载、使用、结果分析方法这篇文章。

1 包的配置与导入

首先，我们可以先到geodetector包在R语言中的官方网站，大致了解一下该包的简要介绍、开发团队、其他依赖包等基本信息；如下图所示。

随后，我们开始geodetector包的下载与安装。输入如下所示的代码，即可开始包的下载与安装过程。

install.packages("geodetector")

输入代码后，按下回车键，运行代码；如下图所示。

随后，将自动下载并配置geodetector包；此外，在安装geodetector包时，会自动将其所需依赖的其他包（如果在此之前没有配置过）都一并配置好，非常方便。

接下来，输入如下的代码，将geodetector包导入。

library(geodetector)

此时，在RStudio右下方的“Packages”中，可以看到geodetector包处于选中的状态，表明其已经配置成功，且完成导入。

2 栅格数据读取与预处理

接下来，我们首先依据基于R语言的raster包读取遥感影像中提到的方法，读取栅格数据。因为我们是要基于栅格数据完成地理探测器的分析，因此很显然是需要批量导入多张栅格数据的。

读取栅格数据完毕后，我们通过如下代码，基于getValues()函数，从原本的RasterStack格式的数据中，将栅格数据的像元数值提取出来；随后，基于View()函数显示出这一变量。

tif_file_all_matrix <- getValues(tif_file_all)
View(tif_file_all_matrix)

运行上述代码，将在RStudio的左上方看到变量tif_file_all_matrix的数据情况，如下图所示。可以看到，此时tif_file_all_matrix变量是一个3列、6377265行的矩阵（Matrix）数据；其中，每一列表示每一个图层的数据，每一行则是每一个图层在同一空间位置上各自像元的数值。此外，每一列的名称即为其所对应的图层的名称。

从上图可以看出，每一列数据中都有很多无效值（NA值），即原本栅格图像中的无效值（NoData值）；由于在后期的地理探测器分析过程中，出现无效值会影响我们分析的结果，因此我们需要通过na.omit()函数将无效值去除。na.omit()是一个非常方便的函数，其可以将Matrix数据中存在NA值的行直接去除（只要这一行中存在至少一个NA，就将这一行去除）。

tif_matrix = na.omit(tif_file_all_matrix)
View(tif_matrix)

随后，我们再看得到的新变量，可以看到存在NA值的行都不复存在了；如下图所示。

接下来，由于geodetector包实现地理探测器操作时，需要保证输入数据为数据框（Data Frames）格式，因此我们需要将Matrix转为Data Frames；通过as.data.frame()函数即可实现这样的转换。

tif_frame <- as.data.frame(tif_matrix)
View(tif_frame)

运行上述代码，可以看到已经获取到Data Frames格式的变量tif_frame了；当然，从外观上看，其和Matrix格式的变量tif_matrix其实长得是一样的。

完成上述数据预处理操作，我们即可开始地理探测器操作。需要注意的是，本文开头也提到了，基于geodetector包实现地理探测器操作时，如果输入的自变量数据是连续数据，我们需要手动将连续数据转为类别数据。这一步骤可以通过ArcGIS的重分类等工具来实现，这里就不再赘述。

3 地理探测器分析

完成上述数据预处理操作，我们即可开始地理探测器的各项具体操作。需要注意的是，本文主要对分析的具体方法加以介绍；至于分析结果的详细研读方法，大家参考文章地理探测器Geodetector下载、使用、结果分析方法即可，我们这里只做简单的介绍。

3.1 分异及因子探测

首先，我们进行分异及因子探测。在geodetector包中，我们可以基于factor_detector()函数实现这一操作。其中，"A_LCCS0"是本文中的因变量，"DEM_Reclass"与"F_LCS0"则是本文中的自变量；tif_frame则是Data Frames格式变量的名称。

在这里需要注意，如果大家只需要分析一个自变量与因变量的影响关系，用下方第一句代码所示的格式即可；如果需要分析多个自变量与因变量的影响关系，则需要用下方第二句代码所示的格式，将多个自变量的名称通过c()函数，组成一个向量（Vector）格式的变量即可。

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

我们首先以上述第一句代码为例来运行，运行后稍等片刻（具体时长与数据量的大小有关），将会得到如下所示的分异及因子探测结果。

其次，再运行上述第二句代码，得到如下所示的结果。

可以看到，factor_detector()函数将会给出每一个自变量对于因变量的q值与p值。

3.2 交互作用探测

接下来，我们执行交互作用探测；这一操作通过interaction_detector()函数来执行即可。由于交互作用探测是需要对多个不同的自变量加以组合，所以很显然这一操作在只有一个自变量的情况下是没有办法执行的；因此我们需要用前述第二种代码格式，即通过c()函数，将多个自变量的名称组成一个向量（Vector）格式的变量后加以执行。

interaction_detector("DEM_Reclass", c("F_LCS0", "K_NDVI"), tif_frame)

运行上述代码，稍等片刻后将出现如下所示的结果。

可以看到，interaction_detector()函数将会给出每一种自变量组合方式对应的q值。但是这里有一个问题——上述结果只能看到不同组合对应的q值变化，但是似乎看不出这种组合方式到底属于非线性减弱、单因子非线性减弱、双因子增强、独立、非线性增强中的哪一种情况。

3.3 风险区探测

接下来，我们执行风险区探测；这一操作通过risk_detector()函数来实现即可，同样是具有一个自变量和多个自变量的情况。我们这里就直接以多个自变量的情况来展示代码与结果了。

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

运行上述代码，稍等片刻后将出现如下所示的结果。

可以看到，risk_detector()函数首先将会给出每一种自变量的不同分级中，对应的因变量平均值——这里自变量的分级指的就是重分类后其的每一个分类；其次，其将给出每一种自变量的分级与分级对应的平均值之间，是否具有显著性差异。

3.4 生态探测

接下来，我们执行生态探测；这一操作通过ecological_detecto()函数来实现即可。由于生态探测是需要判断多个不同的自变量中，两两之间是否具有显著差异，所以很显然这一操作同样在只有一个自变量的情况下是没有办法执行的；因此我们需要用前述第二种代码格式，即通过c()函数，将多个自变量的名称组成一个向量（Vector）格式的变量后加以执行。

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

运行上述代码，稍等片刻后将出现如下所示的结果。

至此，我们就完成了基于R语言中的geodetector包，基于多张栅格图像数据，实现地理探测器（Geodetector）操作的完整流程。