R语言求取大量遥感影像的平均值、标准差：raster库
<p> 本文介绍基于<strong>R</strong>语言中的<code>raster</code>包，批量读取<strong>多张</strong>栅格图像，对多个栅格图像计算<strong>平均值</strong>、<strong>标准差</strong>，并将所得新的栅格结果图像保存的方法。</p><p> 在文章基于R语言的raster包读取遥感影像中，我们介绍了基于<strong>R</strong>语言<code>raster</code>包，对单张或多张栅格图像加以<strong>平均值</strong>、<strong>标准差</strong>计算的方法；但这一篇文章中的<strong>标准差</strong>计算方法仅仅可以对<strong>一张栅格图像</strong>的全部像元加以计算，即标准差计算结果是一个具体的<strong>数值</strong>，而不是一景结果<strong>影像</strong>；无法对多张、多时相的栅格图像进行计算。本文就介绍另一种方法，可以对<strong>多个时相</strong>的大量栅格影像加以逐像元平均值、标准差的计算，从而使得最终的结果是一景表示<strong>各个像元在全部时相的图像中</strong>的平均值或标准差的<strong>图像</strong>。</p>
<p> 首先，我们按照文章基于R语言的raster包读取遥感影像中提到的方法，配置、加载<code>raster</code>包，并通过<code>stack()</code>函数读取同一文件夹下的全部栅格图像，具体代码如下所示。其中，代码的含义我们在上述这一篇文章中已经加以介绍，这里就不再赘述。</p>
<pre><code class="language-r">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)
</code></pre>
<p> 运行上述代码，可以看到已经得到了<code>RasterStack</code>格式的结果数据，如下图所示。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202406/3080295-20240619131457357-1045662349.png"></p>
<p> 接下来，我们通过<code>calc()</code>函数，对多时相栅格遥感影像数据加以计算；其中，其第一个参数<code>tif_file_all</code>就是需要加以计算的多个栅格图像，而第二个参数<code>fun = sd</code>表示我们需要计算标准差；如果我们需要计算平均值，那么就将第二个参数修改为<code>fun = mean</code>即可，我们这里就以标准差为例介绍后续的操作。当然，前述提到的文章基于R语言的raster包读取遥感影像中的方法也是可以对多个栅格图像计算平均值的。</p>
<pre><code class="language-r">tif_sd <- calc(tif_file_all, fun = sd)
plot(tif_sd)
</code></pre>
<p> 此外，上述代码在<code>calc()</code>函数运行时，若某一空间位置上的<strong>像元</strong>在<strong>多张栅格遥感影像</strong>中，存在<strong>至少一个</strong>无效值（<strong>NoData</strong>值），则这一像元在最终的<strong>结果图像</strong>中同样为无效值；若希望忽略无效值的这一影响，可以将上述第一句代码修改为如下格式。其中，<code>na.rm = TRUE</code>就表示若某一景栅格遥感影像中某像元为无效值，则忽略这一景影像中的这一个像元。</p>
<pre><code class="language-r">tif_sd <- calc(tif_file_all, fun = sd, na.rm = TRUE)
</code></pre>
<p> 运行<code>calc()</code>函数后，我们可以通过<code>plot()</code>函数将结果图像绘制出来，如下图所示。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202406/3080295-20240619131451948-384181834.png"></p>
<p> 上图即为多个栅格图像的像元数值时间序列依次计算标准差所得的结果。</p>
<p> 此外，由于我这里的<strong>栅格像元数据</strong>与<strong>实际表达的数值</strong>之间有一个缩放系数<code>0.01</code>，因此通过下述代码将其像元值恢复为实际含义的数值。</p>
<pre><code class="language-r">tif_sd_new <- tif_sd / 100
plot(tif_sd_new)
</code></pre>
<p> 随后，重新绘制结果图；确认无误后，即可依据<code>writeRaster()</code>函数，通过如下代码保存我们刚刚得到的标准差结果栅格图像。</p>
<pre><code class="language-r">rf <- writeRaster(tif_sd_new, filename = r"(E:\02_Project\01_Chlorophyll\LCC_SC_2020\SD\LCC_SD.tif)", overwrite = TRUE)
</code></pre>
<p> 运行代码后，如下图所示。其中，<code>writeRaster()</code>函数的第一个参数表示我们将要保存的栅格数据，第二个参数表示保存栅格文件的路径与名称，第三个参数表示，如果第二个参数指定的路径与名称已经有文件存在了，是否直接对其加以覆盖。</p>
<p><img src="https://img2024.cnblogs.com/blog/3080295/202406/3080295-20240619131452045-673509909.png"></p>
<p> 随后，我们即可在指定的路径下找到我们刚刚计算得到的<strong>多个栅格图像的标准差</strong>结果。</p>
<p> 至此，大功告成。</p><br><br>
来源：https://www.cnblogs.com/fkxxgis/p/18256051
頁:
[1]