（Python）使用Gdal进行批量的多光谱影像波段合成

（Python）使用Gdal进行批量的多光谱影像波段合成

摘要

项目中经常遇到批量多光谱合成的任务需求，数量不多时，可以利用ENVI、ARCGIS等软件进行手工操作，但是当遇到数据量大且数据名称类似的任务时，很容易陷入机械劳动，耗时且易出错。

本文提供了个人学习科研进行中，撰写的批量多光谱合成代码。其可以按照合成对象的命名特点搜索文件中的待合成波段的影像列表，并使用Gdal完成简单的波段合成，最后使用Envi检验合成结果。

方法

代码介绍

代码主要分为两个函数：

1. 自动获待要合成的影像名称和地址。
2. 进行合成图像的创建与多波段的合成。

完整代码

#载入相关库
import os
import re
from osgeo import gdal
import numpy as np

'''
功能：从图中找到要波段融合的影像，然后合成为一个波段（可做批处理）
函数构成：
        1.自动获取待合成的图像名称和地址。
        2.进行图像的创建与多波段的合成。
作者：卫少东
单位：WHU
时间：20211109
'''

def getImageList(image_dir,out_dir,image):
    """
    获取相关待合成的多光谱影像地址列表
    :param image_dir: 影像文件所在文件夹
    :param out_dir: 输出文件夹
    :param image: 待寻找文件名称的正则项
    :return: 创建好的文件列表地址
    """
    image_name = re.compile(image)
    dirPath = image_dir#影像们所在的总目录，也即是遍历进行的目录
    str_size = len(image_name_list[0])
    result_file = out_dir+'\\'+image[0:25]+image[str_size-11:str_size-1]+ r'.txt'#这里按照自己影像的命名规则改变
    fresult_file = open(result_file, 'w')  # 创建影像列表文件
    for root, dirs, files in os.walk(dirPath):
        if files:
            for file in files:
                mo = image_name.search(file)
                if mo:
                    print((root + '\\' + file))
                    fresult_file.write((root + '\\' + file + '\n'))
    print("image list created.")
    return result_file

def SynthesisBands(dstlist,outfile_dir):
    """
    将多光谱波段合成一个tif
    :param dstlist: 输入待合成文件的列表
    :param outfile_dir: 影像的输出文件夹
    """
    image_dst_list = np.loadtxt(dst_list, 'str')
    dataset_init = gdal.Open(image_dst_list[0])
    #创建待输出的图
    gtiff_driver = gdal.GetDriverByName('GTiff')
    dst_name = outfile_dir+ dstlist[dstlist.rfind('\\'):dstlist.rfind('x')] + 'if'#从后往前找第一个"\"注意要用双\\表示
    out_ds = gtiff_driver.Create(dst_name, dataset_init.RasterXSize, dataset_init.RasterYSize, 4)#本文合成的对象是4个波段，按照自己的需要改变
    out_ds.SetProjection(dataset_init.GetProjection())
    out_ds.SetGeoTransform(dataset_init.GetGeoTransform())#获得原始波段的地理信息
    #往图中填各波段
    for i in range(len(image_dst_list)):
        dataset = gdal.Open(image_dst_list[i])
        band_temp = dataset.GetRasterBand(1)
        out_ds.GetRasterBand(1+i).WriteArray(band_temp.ReadAsArray())#注意band从1开始，所以要加一
    del out_ds
    print("the bands systhesised to one tif.")


bands_name = [6,8,14,28] #共四个波段
image_name_list = [r'HHZ3_20201114005237_0007_L1B_B\d\d_CMOS2_ortho',
                   r'HHZ3_20201114005237_0008_L1B_B\d\d_CMOS2_ortho',
                   r'HHZ3_20201114005237_0009_L1B_B\d\d_CMOS2_ortho',
                   r'HHZ3_20201114005237_0010_L1B_B\d\d_CMOS1_ortho',
                   r'HHZ3_20201114005237_0010_L1B_B\d\d_CMOS3_ortho',
                   r'HHZ3_20201114005237_0011_L1B_B\d\d_CMOS1_ortho',
                   r'HHZ3_20201114005237_0011_L1B_B\d\d_CMOS3_ortho'] #这里是正确搜索影像文件名称的正则表达式，需要按照自己的改变
image_dir ="D:\欧比特\欧比特ortho"
out_dir = r'D:\欧比特\欧比特ortho\res'
for name in image_name_list:
    dst_list = getImageList(image_dir, out_dir, name)
    SynthesisBands(dst_list, out_dir)
print("Patch synthesis done.")

实验结果

代码运行结果（部分）

D:\欧比特\欧比特ortho\ortho6\HHZ3_20201114005237_0011_L1B_B06_CMOS3_ortho.tif
D:\欧比特\欧比特ortho\ortho8\HHZ3_20201114005237_0011_L1B_B08_CMOS3_ortho.tif
D:\欧比特\欧比特ortho\ortho14\HHZ3_20201114005237_0011_L1B_B14_CMOS3_ortho.tif
D:\欧比特\欧比特ortho\ortho28\HHZ3_20201114005237_0011_L1B_B28_CMOS3_ortho.tif
image list created.
the bands systhesised to one tif.
Patch synthesis done.

多光谱合成结果

耗时

在商务本ThinkpadX1上，合成一组四张6000*6000左右大小的影像张耗时约3秒。

如需进一步交流，个人联系方式：WHUwsd1995（weixin）