Python遥感开发之批量镶嵌

前言：主要介绍了遥感数据的镶嵌，其中包括使用ArcGis如何完成镶嵌，如何使用Arcpy和GDAL完成镶嵌。

---

1.ArcGis镶嵌

是ArcGis中的非常重要的一个功能，具体Arcgis中的使用可以参考以下博客。
《Arcgis使用教程（五）ARCGIS空间数据处理之影像镶嵌（拼接）与裁剪_arcgis镶嵌工具在哪里-CSDN博客》
简单理解就是由多个TIF合成一个TIF，再合成一个TIF的过程中，可以选择最大值法、中值法、均值法、最小值法等。至于为什么要进行合成，多数是因为单幅影像缺失比较大，需要多幅影像合成才可以，如下图所示！

2.Arcpy实现镶嵌

2.1 Arcpy实现单个镶嵌

其中32_BIT_FLOAT表示浮点型，16_BIT_UNSIGNED表示整型，MAXIMUM表示最大值合成，MEAN表示均值合成。

# encoding:utf-8
import os
from arcpy.management import MosaicToNewRaster

if __name__ == "__main__":

    file = unicode(r"E:\AAWORK\work\b1\20220816", "utf-8")
    out = unicode(r"E:\AAWORK\work\202208合成\bb16arcpy", "utf-8")
    os.chdir(file)
    names = os.listdir(file)
    list1 = []
    for name in names:
        list1.append(file+"\\"+name)
    print list1
    # mosaic = MosaicToNewRaster(list1, out, "202208_b20.tif", pixel_type="32_BIT_FLOAT",number_of_bands=1, mosaic_method="MAXIMUM")
    # mosaic = MosaicToNewRaster(list1, out, "20220816_b4.tif", pixel_type="16_BIT_UNSIGNED",number_of_bands=1, mosaic_method="MAXIMUM")
    mosaic = MosaicToNewRaster(list1, out, "20220816_mean_b1.tif", pixel_type="16_BIT_UNSIGNED",number_of_bands=1, mosaic_method="MEAN")

2.2 Arcpy实现批量镶嵌

在2.1的基础上增加了一个for遍历实现批量处理

# encoding:utf-8
import os
from arcpy.management import MosaicToNewRaster

if __name__ == "__main__":

    name_time = '202108'
    list_band_name = ["b1","b2","b3","b4","b6","b7","b19"]
    out = r"E:\AAWORK\work\{0}合成\bb31mean".format(name_time)
    out = unicode(out, "utf-8")
    for band_name in list_band_name:
        file = r"E:\AAWORK\work\{0}\{1}".format(band_name,name_time)
        file = unicode(file, "utf-8")
        os.chdir(file)
        names = os.listdir(file)
        list1 = []
        for name in names:
            list1.append(file+"\\"+name)
        file_name = "{0}_mean_{1}.tif".format(name_time,band_name)
        mosaic = MosaicToNewRaster(list1, out, file_name, pixel_type="16_BIT_UNSIGNED",number_of_bands=1, mosaic_method="MEAN")
        # mosaic = MosaicToNewRaster(list1, out, file_name, pixel_type="32_BIT_FLOAT",number_of_bands=1, mosaic_method="MEAN")
        print file_name

3.GDAL实现镶嵌

借助GDAL中的ReprojectImage方法实现镶嵌

import os
from osgeo import gdal
from osgeo import gdalconst

def get_data_list(file_path, out=""):
    list1 = []  # 文件的完整路径
    if os.path.isdir(file_path):
        fileList = os.listdir(file_path)
        if out != "":
            for f in fileList:
                out_data = out + "\\" + f
                out_data = out_data.replace(".HDF", "_ndvi.tif")
                list1.append(out_data)
        else:
            for f in fileList:
                pre_data = file_path + '\\' + f  # 文件的完整路径
                list1.append(pre_data)
        return list1


if __name__ == '__main__':
    inputfile_path = r"E:\AAWORK\work\b19\202208"
    output_tif = r"E:\AAWORK\work\202208合成\bb16\20220801-16_b19.tif"
    input_tiffs = get_data_list(inputfile_path)
    # 打开第一个TIFF文件以获取参考信息
    ref_ds = gdal.Open(input_tiffs[0], gdalconst.GA_ReadOnly)
    driver = gdal.GetDriverByName("GTiff")

    # 创建输出TIFF文件
    output_ds = driver.Create(output_tif, ref_ds.RasterXSize, ref_ds.RasterYSize, 1, gdalconst.GDT_Int32)
    output_ds.SetProjection(ref_ds.GetProjection())
    output_ds.SetGeoTransform(ref_ds.GetGeoTransform())

    # 逐一读取和镶嵌TIFF文件
    for input_tiff in input_tiffs:
        input_ds = gdal.Open(input_tiff, gdalconst.GA_ReadOnly)
        band = input_ds.GetRasterBand(1)
        # 镶嵌到输出文件并采用最大值法
        gdal.ReprojectImage(input_ds, output_ds, None, None, gdalconst.GRA_Max)
        # 释放内存
        input_ds = None
    # 关闭输出文件
    output_ds = None
    print("TIFF文件已镶嵌到 %s，采用最大值法。" % output_tif)