R语言批量计算大量栅格图像平均值、标准差

  本文介绍基于R语言中的raster包，批量读取多张栅格图像，对多个栅格图像计算平均值、标准差，并将所得新的栅格结果图像保存的方法。

  在文章R语言读取栅格遥感影像数据的方法（https://www.jianshu.com/p/675f2edb16c3）中，我们介绍了基于R语言raster包，对单张或多张栅格图像加以平均值、标准差计算的方法；但这一篇文章中的标准差计算方法仅仅可以对一张栅格图像的全部像元加以计算，即标准差计算结果是一个具体的数值，而不是一景结果影像；无法对多张、多时相的栅格图像进行计算。本文就介绍另一种方法，可以对多个时相的大量栅格影像加以逐像元平均值、标准差的计算，从而使得最终的结果是一景表示各个像元在全部时相的图像中的平均值或标准差的图像。

  首先，我们按照文章R语言读取栅格遥感影像数据的方法（https://www.jianshu.com/p/675f2edb16c3）中提到的方法，配置、加载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语言读取栅格遥感影像数据的方法（https://www.jianshu.com/p/675f2edb16c3）中的方法也是可以对多个栅格图像计算平均值的。

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()函数的第一个参数表示我们将要保存的栅格数据，第二个参数表示保存栅格文件的路径与名称，第三个参数表示，如果第二个参数指定的路径与名称已经有文件存在了，是否直接对其加以覆盖。

  随后，我们即可在指定的路径下找到我们刚刚计算得到的多个栅格图像的标准差结果。