GDAL多波段合成工具

将一组坐标系相同、行列数相同的单波段tif影像合并为多波段tif影像

先上代码

#!D:\ProgramData\Anaconda3\python.exe
# 指定python指示器

import argparse
import os
import re
from osgeo import gdal
from osgeo import osr
from osgeo import gdal_array
from tqdm import tqdm

# 设置参数
parser = argparse.ArgumentParser(
    '多波段合成工具：将多个单波段tif文件合并为一个多波段tif文件，待合并的tif文件必须有相同的空间参考和行列数')
parser.add_argument('--inputpath', type=str,
                    required=True, help='包含待合并tif文件的目录')
parser.add_argument('--filter', type=str, default=None,
                    help='过滤条件，应当为一个正则表达式，对输入目录中的文件进行筛选')
parser.add_argument('--outpath', type=str, default='./',
                    help='输出文件的目录，默认为当前目录')
parser.add_argument('--outname', type=str, required=True, help='输出文件名')
args = parser.parse_args()

# 获取参数
inputpath = args.inputpath
filterstr = args.filter
outpath = args.outpath
outname = args.outname

# 筛选文件
files = os.listdir(inputpath)
# 只能处理tif格式的文件
files = list(filter(lambda file: file.endswith('.tif'), files))
if len(files) == 0:
    print('只能处理tif文件!')
    exit()
# 基于传入的正则表达式筛选符合条件的文件
files = list(filter(lambda file: re.search(filterstr, file), files))
if len(files) == 0:
    print('没有满足筛选条件的待处理文件!')
    exit()
# 筛选带处理影像的空间参考和行列号
gdal.AllRegister()
srs = None
gt=None
row = -9999
col = -9999
dtypes = []
for file in files:
    ds = gdal.Open(os.path.join(inputpath, file))
    band = gdal.Dataset.GetRasterBand(ds, 1)
    dtypes.append(band.DataType)
    if row == -9999 and col == -9999:
        col = band.XSize
        row = band.YSize
        srs = gdal.Dataset.GetSpatialRef(ds)
        gt=gdal.Dataset.GetGeoTransform(ds)
    else:
        if row != band.YSize or col != band.XSize:
            print('待合并的影像必须有相同的行列号!')
            exit()
        elif not(osr.SpatialReference.IsSame(srs, gdal.Dataset.GetSpatialRef(ds))):
            print('待合并的影像必须有相同的空间参考!')
            exit()
        else:
            pass

driver = gdal.GetDriverByName('GTiff')
N = len(files)
# 设置输出文件的数据类型
dMax = max(dtypes)
# 估算输出文件的大小,基于此决定是否创建BIGTIFF
GDT_dict = dict([
    (1,8),
    (2,16),
    (3,16),
    (4,32),
    (5,32),
    (6,32),
    (7,64),
    (8,16),
    (9,32),
    (10,32),
    (11,64)])
outSize=GDT_dict[dMax]*col*row/8/1024/1024/1024*N
bigTiffFlag='YES' if outSize>=4 else 'NO'
outds = gdal.Driver.Create(driver, os.path.join(
    outpath, outname), col, row, N, dMax, ["COMPRESS=LZW", f"BIGTIFF={bigTiffFlag}"])
gdal.Dataset.SetSpatialRef(outds,srs)
gdal.Dataset.SetGeoTransform(outds,gt)

# 开始合并文件
print('待合并的文件有:\n{}'.format('\n'.join(files)))
pbar = tqdm(files, colour='blue')
for file in pbar:
    pbar.set_description('正在处理{}'.format(file))
    filearr = gdal_array.LoadFile(os.path.join(inputpath, file))
    i=files.index(file)
    outband=gdal.Dataset.GetRasterBand(outds,i+1)
    gdal.Band.WriteArray(outband,filearr)


print('finished!')

示例

bands_composite.py --inputpath=../months_lzw --filter=odiac2020b_1km_excl_intl_00[0-9]{2}.tif --outpath=./ --outname=2000.tif

结果输出

查看输出的多波段tif影像