R语言将多景遥感影像拼接在一起的方法
本文介绍基于R语言中的raster
包,遍历文件夹,读取文件夹下的大量栅格遥感影像,并逐一对每一景栅格图像加以拼接、融合,使得全部栅格遥感影像拼接为完整的一景图像的方法。
其中,本文是用R语言来进行操作的;如果希望基于Python语言实现类似的批量拼接、镶嵌操作,大家可以参考Python arcpy创建栅格、批量拼接栅格与Python ArcPy批量拼接长时间序列栅格图像这两篇文章。
首先,来看一下本文所需实现的需求。如下图所示,现有一个文件夹,其中含有大量栅格遥感影像;这些遥感影像均为同一成像时间、不同空间范围的遥感影像。我们希望做到的,就是对这些遥感影像加以拼接,最终的结果图像就是一景将这里各个图像拼接后的大图像。
明确了需求,我们即可开始代码的撰写。本文所用到的代码如下所示。
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包读取遥感影像。接下来,我们通过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库这篇文章即可。
随后,运行上述代码,我们就可以获得将指定文件夹内全部栅格遥感影像加以拼接(执行代码中的第一部分)或者融合(执行代码中的第二部分)的结果了。
至此,大功告成。