本文介紹基於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）操作的完整流程。