R语言批量拼接遥感影像

2023-06-26 14:47:19 浏览数 (1)

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

  其中,本文是用R语言来进行操作的;如果希望基于Python语言实现类似的批量拼接、镶嵌操作,大家可以参考Python中arcpy栅格创建与多景遥感影像批量拼接与ArcPy依据成像时间分别批量拼接不同时相的遥感影像这两篇文章。

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

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

代码语言:javascript复制
library(raster)
tif_file_name <- list.files(path = r"(E:2_Project1_ChlorophyllSelectResult)", 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:2_Project1_ChlorophyllSelectNewClipLCC_SC_3.tif)", overwrite = TRUE)

  首先,需要通过library(raster)代码,导入本文所需的R语言raster包;关于这一包的配置,大家可以参考R语言raster包读取栅格遥感影像。接下来,我们通过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个参数分别设定为maxTRUE,二者分别是指重叠区域的像元以2景遥感影像中的最大值像元为准,以及在计算时不考虑NoData值的影响。

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

  最后,这里末尾的一句代码,就是将结果图像通过writeRaster()函数加以保存;这句代码的解释大家同样参考R语言计算大量栅格图像平均值、标准差这篇文章即可。

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

  至此,大功告成。

0 人点赞