R语言GD包基于栅格图像实现地理探测器与连续参数的自动离散化

  本文介绍基于R语言中的GD包，依据栅格影像数据，实现自变量最优离散化方法选取与执行，并进行地理探测器（Geodetector）操作的方法。

  首先，在R语言中进行地理探测器操作，可通过geodetector包、GD包等2个包实现。其中，geodetector包是地理探测器模型的原作者团队开发的，其需要保证输入的自变量数据已经全部为类别数据；其具体操作方法大家可以参考R语言geodetector包基于栅格图像实现地理探测器操作（https://blog.csdn.net/zhebushibiaoshifu/article/details/128484786）。而GD包则是另一位学者开发的，其可自动实现自变量数据的最优离散化方法选取与执行；本文介绍的就是基于GD包实现地理探测器的具体操作。此外，如果希望基于Excel实现地理探测器，大家可以参考Geodetector软件下载、地理探测器的应用实践与结果解读（https://blog.csdn.net/zhebushibiaoshifu/article/details/124526189）。

1 包的安装与导入

  首先，我们可以先到GD包在R语言中的官方网站（https://cran.r-project.org/web/packages/GD/index.html），大致了解一下该包的简要介绍、开发团队等基本信息。

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

install.packages("GD")

  输入代码后，按下回车键，运行代码；如下图所示。在安装GD包时，会自动将其所需依赖的其他包（如果在此之前没有配置过）都一并配置好，非常方便。

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

library("GD")

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

2 数据读取与预处理

  接下来，我们需要读取栅格图像数据，并将其转为GD包可以识别的数据框（Data Frames）格式。

  其中，读取栅格数据的方法，大家参考R语言raster包批量读取单一或大量栅格图像（https://blog.csdn.net/zhebushibiaoshifu/article/details/128485386）即可；关于数据格式的转换，大家参考R语言geodetector包基于栅格图像实现地理探测器操作（https://blog.csdn.net/zhebushibiaoshifu/article/details/128484786）即可。这一部分的内容本文就不再赘述。

3 地理探测器执行

  接下来，我们就可以开始地理探测器的具体分析；强烈建议大家基于GD包中的gdm()函数，实现一步到位的地理探测器分析操作。

  首先，如果大家输入数据中的自变量数据具有连续变量，需要将其转换为类别变量；gdm()函数可以实现连续变量的离散化方式寻优与自动执行。其中，我们可以选择的离散化方式包括相等间隔法、自然间断点法、分位数分类法、几何间隔法与标准差法等5种不同的方法，分别对应以下第一句代码中的"equal"、"natural"、"quantile"、"geometric"与"sd"等5个选项。此外，我们还可以依据数据的特征，对自变量离散化的类别数量加以限定，具体代码如下所示。

discmethod <- c("equal", "natural", "quantile", "geometric", "sd")
discitv <- c(4:10)

  其中，上述第一句代码表示，我们后续将从相等间隔法、自然间断点法、分位数分类法、几何间隔法与标准差法等5种不同的方法中，找到每一个连续变量对应的最优离散化方法；第二句代码则表示，在后续寻找最优离散化方法的同时，还需要对每一个变量的分类数量加以寻优——c(4:10)就表示我们分别将每一个连续变量分为4类、5类、6类，以此类推，一直到10类，从其中找到最优结果对应的类别数量。

  接下来，我们即可调用gdm()函数，执行地理探测器分析的具体操作；其中，my_gd为保存地理探测器结果的变量；函数的第一个参数，表示因变量与自变量的关系，~前的变量即为因变量，~后的变量即为自变量，多个自变量之间通过+相连接；第二个参数表示自变量中的连续变量，程序将自动对这些连续变量加以离散化方法寻优与执行；第三个参数表示存储自变量与因变量数据的数据框（Data Frames）格式的变量；最后两个变量，即为前面我们选择的离散化方法与类别数量。

my_gd <- gdm(A_LCCS0 ~ C_SlopeS0 + D_AspectS0 + DEM_Reclass + F_LCS0,
                        continuous_variable = c("C_SlopeS0", "D_AspectS0"),
                        data = tif_frame,
                        discmethod = discmethod,
                        discitv = discitv)

  这里需要注意，如果大家不是通过脚本运行的R语言，而是每次写一句代码然后按下回车键运行一下，那么上述代码中的换行就需要通过同时按下Shift键与回车键实现。输入上述代码后，如下图所示。

  随后，即可运行代码。稍等片刻（具体时长与数据量有关），即可得到地理探测器的结果my_gd。这一变量的具体结构、内容如下图所示。

  我们可以输入如下的代码，将变量my_gd打印出来。

my_gd

  所得结果如下图所示。

  可以看到，my_gd变量包含了每一个连续变量在离散化后，对应的最优离散化方法与类别数量，以及地理探测器的各个分析结果。具体结果的含义与研读方法，大家参考Geodetector软件下载、地理探测器的应用实践与结果解读（https://blog.csdn.net/zhebushibiaoshifu/article/details/124526189）、R语言geodetector包基于栅格图像实现地理探测器操作（https://blog.csdn.net/zhebushibiaoshifu/article/details/128484786）这两篇文章即可，这里就不再赘述。

  此外，我们可以通过如下的代码，将上述结果加以可视化。

plot(my_gd)

  运行上述代码，结果如下图所示。

  此时，在RStudio软件的右下方“Plots”中，即可看到可视化结果，如下图所示。其中，我们可以通过下图中红色方框内的箭头，实现不同图片的切换显示。

  上述结果包含7张图像，其分别与上上图中的7项输出内容对应——第一张图是最优离散化方法的选取过程，第二张图则是所选出的最优离散化方法对应的分类情况；后5张图就是地理探测器的分析结果图，即上上图中最后5个plot分别对应的结果。

  至此，我们就完成了基于R语言中的GD包，依据多张栅格图像数据，实现类别变量的自动离散化，并进行地理探测器（Geodetector）操作的完整流程。

欢迎关注：疯狂学习GIS