渐难医 發表於 2024-3-18 10:28:00

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

<p>  本文介绍基于<strong>R</strong>语言中的<code>geodetector</code>包,依据多张<strong>栅格图像</strong>数据,实现<strong>地理探测器</strong>(<strong>Geodetector</strong>)操作的详细方法。</p>
<p>  需要说明的是,在<strong>R</strong>语言中进行<strong>地理探测器</strong>操作,可以分别通过<code>geodetector</code>包、<code>GD</code>包等<code>2</code>个包实现。其中,<code>geodetector</code>包是<strong>地理探测器模型</strong>的原作者团队早先开发的,其需要保证输入的<strong>自变量数据</strong>已经全部为<strong>类别数据</strong>;而<code>GD</code>包则是另外一位学者开发的,其可以自动实现自变量数据的<strong>最优离散化方法</strong>选取与执行——即我们可以直接把自变量带入这一包中,无需额外进行数据的离散化。本文介绍的是基于前者,即<code>geodetector</code>包实现地理探测器的具体操作;基于后者的方法,我们将在后期的博客中介绍。此外,如果希望基于<strong>Excel</strong>实现<strong>地理探测器</strong>,大家可以参考地理探测器Geodetector下载、使用、结果分析方法这篇文章。</p>
<h1 id="1-包的配置与导入">1 包的配置与导入</h1>
<p>  首先,我们可以先到<code>geodetector</code>包在<strong>R</strong>语言中的官方网站,大致了解一下该包的简要介绍、开发团队、其他依赖包等基本信息;如下图所示。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202403/3080295-20240318102721908-1949109217.png"></p>
<p>  随后,我们开始<code>geodetector</code>包的下载与安装。输入如下所示的代码,即可开始包的下载与安装过程。</p>
<pre><code class="language-r">install.packages("geodetector")
</code></pre>
<p>  输入代码后,按下<code>回车</code>键,运行代码;如下图所示。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202403/3080295-20240318102713960-1134465858.png"></p>
<p>  随后,将自动下载并配置<code>geodetector</code>包;此外,在安装<code>geodetector</code>包时,会自动将其所需依赖的其他包(如果在此之前没有配置过)都一并配置好,非常方便。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202403/3080295-20240318102713974-1871064221.png"></p>
<p>  接下来,输入如下的代码,将<code>geodetector</code>包导入。</p>
<pre><code class="language-r">library(geodetector)
</code></pre>
<p>  此时,在<strong>RStudio</strong>右下方的“<strong>Packages</strong>”中,可以看到<code>geodetector</code>包处于选中的状态,表明其已经配置成功,且完成导入。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202403/3080295-20240318102714078-1013220263.png"></p>
<h1 id="2-栅格数据读取与预处理">2 栅格数据读取与预处理</h1>
<p>  接下来,我们首先依据基于R语言的raster包读取遥感影像中提到的方法,读取栅格数据。因为我们是要基于栅格数据完成地理探测器的分析,因此很显然是需要批量导入多张栅格数据的。</p>
<p>  读取栅格数据完毕后,我们通过如下代码,基于<code>getValues()</code>函数,从原本的<code>RasterStack</code>格式的数据中,将栅格数据的像元数值提取出来;随后,基于<code>View()</code>函数显示出这一变量。</p>
<pre><code class="language-r">tif_file_all_matrix &lt;- getValues(tif_file_all)
View(tif_file_all_matrix)
</code></pre>
<p>  运行上述代码,将在<strong>RStudio</strong>的左上方看到变量<code>tif_file_all_matrix</code>的数据情况,如下图所示。可以看到,此时<code>tif_file_all_matrix</code>变量是一个<code>3</code>列、<code>6377265</code>行的<strong>矩阵</strong>(<code>Matrix</code>)数据;其中,每一列表示每一个图层的数据,每一行则是每一个图层在同一空间位置上各自像元的数值。此外,每一列的名称即为其所对应的图层的名称。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202403/3080295-20240318102714016-1994711771.png"></p>
<p>  从上图可以看出,每一列数据中都有很多<strong>无效值</strong>(<strong>NA</strong>值),即原本栅格图像中的<strong>无效值</strong>(<strong>NoData</strong>值);由于在后期的地理探测器分析过程中,出现无效值会影响我们分析的结果,因此我们需要通过<code>na.omit()</code>函数将无效值去除。<code>na.omit()</code>是一个非常方便的函数,其可以将<code>Matrix</code>数据中存在<strong>NA</strong>值的行直接去除(只要这一行中存在至少一个<strong>NA</strong>,就将这一行去除)。</p>
<pre><code class="language-r">tif_matrix = na.omit(tif_file_all_matrix)
View(tif_matrix)
</code></pre>
<p>  随后,我们再看得到的新变量,可以看到存在<strong>NA</strong>值的行都不复存在了;如下图所示。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202403/3080295-20240318102713922-970328839.png"></p>
<p>  接下来,由于<code>geodetector</code>包实现地理探测器操作时,需要保证输入数据为<strong>数据框</strong>(<code>Data Frames</code>)格式,因此我们需要将<code>Matrix</code>转为<code>Data Frames</code>;通过<code>as.data.frame()</code>函数即可实现这样的转换。</p>
<pre><code class="language-r">tif_frame &lt;- as.data.frame(tif_matrix)
View(tif_frame)
</code></pre>
<p>  运行上述代码,可以看到已经获取到<code>Data Frames</code>格式的变量<code>tif_frame</code>了;当然,从外观上看,其和<code>Matrix</code>格式的变量<code>tif_matrix</code>其实长得是一样的。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202403/3080295-20240318102714093-1632599221.png"></p>
<p>  完成上述数据预处理操作,我们即可开始地理探测器操作。需要注意的是,本文开头也提到了,基于<code>geodetector</code>包实现地理探测器操作时,如果输入的自变量数据是连续数据,我们需要手动将连续数据转为类别数据。这一步骤可以通过<strong>ArcGIS</strong>的重分类等工具来实现,这里就不再赘述。</p>
<h1 id="3-地理探测器分析">3 地理探测器分析</h1>
<p>  完成上述数据预处理操作,我们即可开始地理探测器的各项具体操作。需要注意的是,本文主要对分析的具体方法加以介绍;至于<strong>分析结果的详细研读方法</strong>,大家参考文章地理探测器Geodetector下载、使用、结果分析方法即可,我们这里只做简单的介绍。</p>
<h2 id="31-分异及因子探测">3.1 分异及因子探测</h2>
<p>  首先,我们进行<strong>分异及因子探测</strong>。在<code>geodetector</code>包中,我们可以基于<code>factor_detector()</code>函数实现这一操作。其中,<code>"A_LCCS0"</code>是本文中的因变量,<code>"DEM_Reclass"</code>与<code>"F_LCS0"</code>则是本文中的自变量;<code>tif_frame</code>则是<code>Data Frames</code>格式变量的名称。</p>
<p>  在这里需要注意,如果大家只需要分析<strong>一个自变量</strong>与因变量的影响关系,用下方第一句代码所示的格式即可;如果需要分析<strong>多个自变量</strong>与因变量的影响关系,则需要用下方第二句代码所示的格式,将多个自变量的名称通过<code>c()</code>函数,组成一个<strong>向量</strong>(<code>Vector</code>)格式的变量即可。</p>
<pre><code class="language-r">factor_detector("A_LCCS0", "F_LCS0", tif_frame)
factor_detector("A_LCCS0", c("DEM_Reclass", "F_LCS0"), tif_frame)
</code></pre>
<p>  我们首先以上述第一句代码为例来运行,运行后稍等片刻(具体时长与数据量的大小有关),将会得到如下所示的<strong>分异及因子探测</strong>结果。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202403/3080295-20240318102714034-1180915730.png"></p>
<p>  其次,再运行上述第二句代码,得到如下所示的结果。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202403/3080295-20240318102714142-91831926.png"></p>
<p>  可以看到,<code>factor_detector()</code>函数将会给出每一个自变量对于因变量的<code>q</code>值与<code>p</code>值。</p>
<h2 id="32-交互作用探测">3.2 交互作用探测</h2>
<p>  接下来,我们执行<strong>交互作用探测</strong>;这一操作通过<code>interaction_detector()</code>函数来执行即可。由于<strong>交互作用探测</strong>是需要对多个不同的自变量加以组合,所以很显然这一操作在只有<strong>一个自变量</strong>的情况下是没有办法执行的;因此我们需要用前述第二种代码格式,即通过<code>c()</code>函数,将多个自变量的名称组成一个<strong>向量</strong>(<code>Vector</code>)格式的变量后加以执行。</p>
<pre><code class="language-r">interaction_detector("DEM_Reclass", c("F_LCS0", "K_NDVI"), tif_frame)
</code></pre>
<p>  运行上述代码,稍等片刻后将出现如下所示的结果。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202403/3080295-20240318102714093-1464262468.png"></p>
<p>  可以看到,<code>interaction_detector()</code>函数将会给出每一种自变量组合方式对应的<code>q</code>值。但是这里有一个问题——上述结果只能看到不同组合对应的<code>q</code>值变化,但是似乎看不出这种组合方式到底属于<strong>非线性减弱</strong>、<strong>单因子非线性减弱</strong>、<strong>双因子增强</strong>、<strong>独立</strong>、<strong>非线性增强</strong>中的哪一种情况。</p>
<h2 id="33-风险区探测">3.3 风险区探测</h2>
<p>  接下来,我们执行<strong>风险区探测</strong>;这一操作通过<code>risk_detector()</code>函数来实现即可,同样是具有<strong>一个自变量</strong>和<strong>多个自变量</strong>的情况。我们这里就直接以多个自变量的情况来展示代码与结果了。</p>
<pre><code class="language-r">risk_detector("A_LCCS0", c("DEM_Reclass", "F_LCS0"), tif_frame)
</code></pre>
<p>  运行上述代码,稍等片刻后将出现如下所示的结果。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202403/3080295-20240318102714110-1740005843.png"></p>
<p>  可以看到,<code>risk_detector()</code>函数首先将会给出每一种自变量的<strong>不同分级</strong>中,对应的因变量平均值——这里自变量的<strong>分级</strong>指的就是重分类后其的<strong>每一个分类</strong>;其次,其将给出每一种自变量的分级与分级<strong>对应的平均值</strong>之间,是否具有显著性差异。</p>
<h2 id="34-生态探测">3.4 生态探测</h2>
<p>  接下来,我们执行<strong>生态探测</strong>;这一操作通过<code>ecological_detecto()</code>函数来实现即可。由于<strong>生态探测</strong>是需要判断多个不同的自变量中,两两之间是否具有显著差异,所以很显然这一操作同样在只有一个自变量的情况下是没有办法执行的;因此我们需要用前述第二种代码格式,即通过<code>c()</code>函数,将多个自变量的名称组成一个<strong>向量</strong>(<code>Vector</code>)格式的变量后加以执行。</p>
<pre><code class="language-r">ecological_detector("A_LCCS0", c("DEM_Reclass", "F_LCS0"), tif_frame)
</code></pre>
<p>  运行上述代码,稍等片刻后将出现如下所示的结果。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202403/3080295-20240318102714079-1198615043.png"></p>
<p>  至此,我们就完成了基于<strong>R</strong>语言中的<code>geodetector</code>包,基于多张<strong>栅格图像</strong>数据,实现<strong>地理探测器</strong>(<strong>Geodetector</strong>)操作的完整流程。</p><br><br>
来源:https://www.cnblogs.com/fkxxgis/p/18079822
頁: [1]
查看完整版本: 地理探测器R语言实现:geodetector