R语言raster包批量拼接、融合大量栅格图像

  本文介绍基于R语言中的raster包，遍历文件夹，读取文件夹下的大量栅格遥感影像，并逐一对每一景栅格图像加以拼接、融合，使得全部栅格遥感影像拼接为完整的一景图像的方法。

  其中，本文是用R语言来进行操作的；如果希望基于Python语言实现类似的批量拼接、镶嵌操作，大家可以参考Python中arcpy栅格创建与遥感影像多景数据批量拼接Mosaic（https://blog.csdn.net/zhebushibiaoshifu/article/details/118901707）与Python中ArcPy实现对不同时相的栅格遥感影像依据其成像时间分别批量拼接（https://blog.csdn.net/zhebushibiaoshifu/article/details/124227237）这两篇文章。

  首先，来看一下本文所需实现的需求。如下图所示，现有一个文件夹，其中含有大量栅格遥感影像；这些遥感影像均为同一成像时间、不同空间范围的遥感影像。我们希望做到的，就是对这些遥感影像加以拼接，最终的结果图像就是一景将这里各个图像拼接后的大图像。

  明确了需求，我们即可开始代码的撰写。本文所用到的代码如下所示。

library(raster)
tif_file_name <- list.files(path = r"(E:\02_Project\01_Chlorophyll\Select\Result)", pattern = ".tif$", full.names = TRUE, ignore.case = TRUE)
tif_file_list <- list()
for (i in 1:length(tif_file_name)){
  tif_file_list[i] <- raster(tif_file_name[i])
}

tif_file_list$fun <- max
tif_file_list$na.rm <- TRUE
tif_mosaic <- do.call(mosaic, tif_file_list)
plot(tif_mosaic)

# tif_merge <- do.call(merge, tif_file_list)

rf <- writeRaster(tif_mosaic, filename = r"(E:\02_Project\01_Chlorophyll\Select\NewClip\LCC_SC_3.tif)", overwrite = TRUE)

  首先，需要通过library(raster)代码，导入本文所需的R语言raster包；关于这一包的配置，大家可以参考R语言raster包批量读取单一或大量栅格图像（https://blog.csdn.net/zhebushibiaoshifu/article/details/128485386）。接下来，我们通过list.files()函数，遍历指定文件夹，从而获取当前文件夹下所包含的全部.tif格式的遥感影像，也就是全部待拼接的遥感影像。

  接下来，我们需要为栅格遥感影像的拼接做准备——也就是for循环内部的内容。此时，tif_file_name变量中存放的是指定文件夹下的全部栅格遥感影像的文件名称，而不是遥感影像文件自身；而接下来我们进行拼接、融合的函数，都需要保证函数参数中的遥感影像是一个栅格对象（Raster* object）类型的变量。因此，我们需要在这个for循环中，通过raster()函数，将每一个遥感影像的文件名（字符串类型）转为栅格对象类型。至于什么是栅格对象类型的变量，我们可以参考下图：其中Formal class RasterLayer即表示这一变量为栅格对象类型的。

  接下来，代码分为2个部分。其中，for循环后的4行代码是第一部分，为栅格拼接的代码；同时为了对比栅格拼接与栅格融合的操作，这里还将栅格融合的代码也一并列出了，也就是注释掉的那一行代码。

  我们首先来看第一部分代码，这里通过mosaic()函数来实现栅格遥感影像的拼接。这一函数原本的参数中，只有2个栅格对象（Raster* object）类型的参数，换句话说就是原本这个函数只能同时拼接2个栅格遥感影像；如果我们有更多的遥感影像，就需要每一次拼接2个栅格图像，不断重复这一操作，直到全部的栅格遥感影像拼接完毕。这样操作无疑是比较麻烦的，因此我们需要借助do.call()函数来实现2个以上栅格的拼接工作——这个do.call()函数可以接受可变数量的参数，例如本文中我们需要对大量栅格遥感影像加以逐一拼接，具体有多少景遥感影像我们自己也不一定确定，且也不关心；因此就结合这一函数，将刚刚已经转为栅格对象（Raster* object）类型的图像所组成的列表tif_file_list作为参数，用do.call()函数来调用mosaic()函数，直到将tif_file_list列表中全部的栅格对象（Raster* object）类型的元素都带入到mosaic()函数运行后，do.call()函数就结束了。

  此外，由于mosaic()函数在运行时，除了两个栅格对象（Raster* object）类型的参数，还有其他的一些辅助参数，比如拼接时重叠区域该如何处理、处理时是否考虑NoData值的影响等；由于我们时通过do.call()函数来调用mosaic()函数，因此这些参数就不太好直接指定了。因此，我们可以通过$运算符，将mosaic()函数所需要的其他参数一并放入tif_file_list中，在后期do.call()函数调用mosaic()函数时，将同时读取这些参数，起到将参数传递到mosaic()函数中的功能。其中，在本文中我们需要指定mosaic()函数的fun参数与na.rm参数，二者分别是指拼接时重叠区域像元值的计算方法，以及计算重叠区域像元值时，是否考虑NoData值的影响；我们将这2个参数分别设定为max与TRUE，二者分别是指重叠区域的像元以2景遥感影像中的最大值像元为准，以及在计算时不考虑NoData值的影响。

  接下来，就是第二部分，即栅格融合的代码；在这里，我们通过merge()函数来实现遥感影像的融合。其实，这里的merge()函数与前述的mosaic()函数功能大致一样，但merge()函数在处理重叠区域时，默认选择位于顶层的遥感影像的像元数值，就没有mosaic()函数中的这么多计算方法选择了。

  最后，这里末尾的一句代码，就是将结果图像通过writeRaster()函数加以保存；这句代码的解释大家同样参考R语言raster包计算多个栅格图像平均值、标准差的方法（https://blog.csdn.net/zhebushibiaoshifu/article/details/128657182）这篇文章即可。

  随后，运行上述代码，我们就可以获得将指定文件夹内全部栅格遥感影像加以拼接（执行代码中的第一部分）或者融合（执行代码中的第二部分）的结果了。

  至此，大功告成。

欢迎关注：疯狂学习GIS