天天看點

R語言raster包周遊檔案夾并計算大量栅格資料平均值

作者:瘋狂學習GIS

  本文介紹基于R語言中的raster包,周遊讀取多個檔案夾下的多張栅格遙感影像,分别批量對每一個檔案夾中的多個栅格圖像計算平均值,并将所得各個結果栅格分别加以儲存的方法。

  其中,本文是用R語言來進行操作的;如果希望基于Python語言實作類似的平均值求取操作,大家可以參考ArcPy對多張長時間序列栅格影像逐像元批量求取像素平均值與Python中whitebox忽略無效值NoData求取多時相栅格圖像像素平均值這兩篇文章。

  首先,來看一下本文所需實作的需求。如下圖所示,現有多個檔案夾,其中每一個檔案夾内部都含有大量的栅格遙感影像。

R語言raster包周遊檔案夾并計算大量栅格資料平均值

  其中,上圖中的每一個檔案夾的命名都是以遙感影像的分幅條帶号為依據的。例如,打開第一個名為47RMN的檔案夾,其中均為條帶号為47RMN(即同一空間範圍)、不同成像時間的遙感影像,如下圖所示;其中,紫色框内的遙感影像檔案名即可看出,這些圖像是同一條帶号、不同時間的遙感影像資料。

R語言raster包周遊檔案夾并計算大量栅格資料平均值

  我們要做的,就是分别對每一個檔案夾中的全部遙感影像計算平均值,進而得到不同條帶号遙感影像的平均值;最終我們将得到多張結果圖像,每一景結果圖像就是這一條帶号、不同成像時間對應的遙感影像的平均值。同時為了友善區分,我們需要将每一景結果圖像檔案的檔案名設定為與條帶号有關的内容。

  明确了需求,我們即可開始代碼的撰寫。本文所用到的代碼如下所示。

library(raster)
result_path <- r"(E:\02_Project\01_Chlorophyll\Select\Result)"
tif_folder <- list.files(path = r"(E:\02_Project\01_Chlorophyll\Select)", pattern = NULL, all.files = FALSE, full.names = TRUE)
for (folder in tif_folder){
  folder_name <- substr(folder, nchar(folder) - 4, nchar(folder))
  tif_file_name <- list.files(path = folder, pattern = ".tif#34;, full.names = TRUE, ignore.case = TRUE)
  tif_file_all <- stack(tif_file_name)
  NAvalue(tif_file_all) <- -10000
  tif_mean <- calc(tif_file_all, fun = mean, na.rm = TRUE)
  tif_mean_new <- tif_mean / 100
  # plot(tif_mean_new)
  result_file_name <- file.path(result_path, paste(folder_name, "_mean.tif", sep = ""))
  rf <- writeRaster(tif_mean_new, filename = result_file_name, overwrite = TRUE)
  cat(folder_name, "is completed!", "\n")
}           

  首先,需要通過library(raster)代碼,導入本文所需的R語言raster包;關于這一包的配置,大家可以參考R語言讀取栅格遙感影像資料的方法。接下來,我們需要指定結果存放的路徑,并将其放入變量result_path中。

  接下來,我們通過list.files()函數,将包含有各個條帶号的小檔案夾的大檔案夾(也就是本文開頭第一張圖所示的檔案夾)加以周遊,将每一個小檔案夾的路徑存入tif_folder。執行上述前3行代碼後,得到的tif_folder結果如下圖所示。

R語言raster包周遊檔案夾并計算大量栅格資料平均值

  可以看到,tif_folder是一個字元串,其中每一個元素都是每一個小檔案夾的路徑。

  接下來的for循環,就是對tif_folder加以周遊,即對每一個小檔案夾進行操作。其中,我們首先通過substr()函數,擷取目前操作的小檔案夾名稱,并将其存放于folder_name中;随後,對目前對應的小檔案夾加以周遊,取出其中的全部遙感影像檔案,并存放于tif_file_name;接下來,就是讀取全部遙感影像,并計算其平均值;這裡具體的代碼解釋大家可以參考文章R語言計算遙感影像的平均值、标準差。此外需要注意的是,由于我這裡每一景遙感影像原本沒有專門設定NoData值,而是用-10000作為其NoData值,是以需要通過NAvalue(tif_file_all) <- -10000這句代碼,将值為-10000的像元作為NoData值的像元,防止後期計算平均值時對結果加以幹擾。

  接下來,我們通過file.path()函數配置一下輸出結果的路徑——其中,結果遙感影像檔案的名稱就可以直接以其所對應的條帶号來設定,并在條帶号後添加一個_mean字尾,表明這個是平均值的結果圖像;但此外,這個僅僅是檔案的名字,還需要将檔案名與路徑拼接在一起,才可以成為完整的儲存路徑,是以需要用到file.path()函數。最後,将結果圖像通過writeRaster()函數加以儲存即可,這句代碼的解釋大家同樣參考R語言計算遙感影像的平均值、标準差這篇文章即可。

  最後,由于我們要處理的檔案夾比較多,是以可以通過cat()函數輸出一下目前代碼的運作進度。

  運作上述代碼,我們将在指定的結果儲存路徑中看到每一個條帶号對應的平均值結果圖像,如下圖所示。

R語言raster包周遊檔案夾并計算大量栅格資料平均值

  至此,大功告成。

歡迎關注:瘋狂學習GIS

R語言raster包周遊檔案夾并計算大量栅格資料平均值

繼續閱讀