本文介紹基於R語言中的raster
包,批次讀取多張柵格影像,對多個柵格影像計算平均值、標準差,並將所得新的柵格結果影像儲存的方法。
在文章基於R語言的raster包讀取遙感影像中,我們介紹了基於R語言raster
包,對單張或多張柵格影像加以平均值、標準差計算的方法;但這一篇文章中的標準差計算方法僅僅可以對一張柵格影像的全部像元加以計算,即標準差計算結果是一個具體的數值,而不是一景結果影像;無法對多張、多時相的柵格影像進行計算。本文就介紹另一種方法,可以對多個時相的大量柵格影像加以逐像元平均值、標準差的計算,從而使得最終的結果是一景表示各個像元在全部時相的影像中的平均值或標準差的影像。
首先,我們按照文章基於R語言的raster包讀取遙感影像中提到的方法,配置、載入raster
包,並透過stack()
函式讀取同一資料夾下的全部柵格影像,具體程式碼如下所示。其中,程式碼的含義我們在上述這一篇文章中已經加以介紹,這裡就不再贅述。
library(raster)
tif_file_path <- list.files(r"(E:\02_Project\01_Chlorophyll\LCC_SC_2020\SD)", pattern = ".tif$", full.names = TRUE, ignore.case = TRUE)
tif_file_all <- stack(tif_file_path)
執行上述程式碼,可以看到已經得到了RasterStack
格式的結果資料,如下圖所示。
接下來,我們透過calc()
函式,對多時相柵格遙感影像資料加以計算;其中,其第一個引數tif_file_all
就是需要加以計算的多個柵格影像,而第二個引數fun = sd
表示我們需要計算標準差;如果我們需要計算平均值,那麼就將第二個引數修改為fun = mean
即可,我們這裡就以標準差為例介紹後續的操作。當然,前述提到的文章基於R語言的raster包讀取遙感影像中的方法也是可以對多個柵格影像計算平均值的。
tif_sd <- calc(tif_file_all, fun = sd)
plot(tif_sd)
此外,上述程式碼在calc()
函式執行時,若某一空間位置上的像元在多張柵格遙感影像中,存在至少一個無效值(NoData值),則這一像元在最終的結果影像中同樣為無效值;若希望忽略無效值的這一影響,可以將上述第一句程式碼修改為如下格式。其中,na.rm = TRUE
就表示若某一景柵格遙感影像中某像元為無效值,則忽略這一景影像中的這一個像元。
tif_sd <- calc(tif_file_all, fun = sd, na.rm = TRUE)
執行calc()
函式後,我們可以透過plot()
函式將結果影像繪製出來,如下圖所示。
上圖即為多個柵格影像的像元數值時間序列依次計算標準差所得的結果。
此外,由於我這裡的柵格像後設資料與實際表達的數值之間有一個縮放係數0.01
,因此透過下述程式碼將其像元值恢復為實際含義的數值。
tif_sd_new <- tif_sd / 100
plot(tif_sd_new)
隨後,重新繪製結果圖;確認無誤後,即可依據writeRaster()
函式,透過如下程式碼儲存我們剛剛得到的標準差結果柵格影像。
rf <- writeRaster(tif_sd_new, filename = r"(E:\02_Project\01_Chlorophyll\LCC_SC_2020\SD\LCC_SD.tif)", overwrite = TRUE)
執行程式碼後,如下圖所示。其中,writeRaster()
函式的第一個參數列示我們將要儲存的柵格資料,第二個參數列示儲存柵格檔案的路徑與名稱,第三個參數列示,如果第二個引數指定的路徑與名稱已經有檔案存在了,是否直接對其加以覆蓋。
隨後,我們即可在指定的路徑下找到我們剛剛計算得到的多個柵格影像的標準差結果。
至此,大功告成。