R语言求取大量遥感影像的平均值、标准差：raster库

raster · 浏览次数 : 0

小编点评

本文介绍了一种基于R语言的raster包批量读取多张栅格图像，并对多个栅格图像计算平均值、标准差的方法。文章首先概述了基于R语言的raster包读取遥感影像的方法，然后详细介绍了如何对多时相的大量栅格影像进行逐像元平均值、标准差的计算，最后展示了如何将计算结果保存为栅格图像。 1. 文章首先介绍了基于R语言的raster包读取遥感影像的方法，包括配置、加载raster包，以及通过stack()函数读取同一文件夹下的全部栅格图像。 - 配置和加载raster包 - 使用stack()函数读取栅格图像 2. 文章接着介绍了如何对多时相栅格遥感影像数据计算平均值和标准差，以及如何处理多张栅格图像中存在的无效值。 - 计算平均值和标准差的方法 - 处理多张栅格图像中的无效值 3. 文章展示了如何通过plot()函数将计算结果图像绘制出来，并将结果图像保存为栅格图像。 - 绘制结果图像 - 保存结果栅格图像 4. 文章最后介绍了如何将栅格像元值恢复为实际含义的数值，并重新绘制结果图，确认无误后，即可依据writeRaster()函数，将计算得到的标准差结果栅格图像保存到指定路径。总的来说，本文提供了一种高效的方法来处理大量栅格图像的平均值和标准差计算，并将结果以栅格图像的形式保存，便于分析和可视化。

正文

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

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

至此，大功告成。