【Pyhton】多年dat数据波段合成为一个dat

我有一个文件夹，内有1961-2020年逐年平均气温栅格数据，格式为dat，dat文件不是以年份命名的，而是一长串名字的最后4位是年份，现在我要对所有年份进行波段合成，合成一个dat，dat每个波段为不同年份的数据，且波段顺序按照年份排列，生成的dat文件能够在ENVI中打开。

import os
from osgeo import gdal

# 设置工作目录和输出文件路径
input_folder = r"D:\itm\平均气温"
output_file = r"D:\itm\平均气温.dat"

# 创建一个列表用于存储每年的栅格数据路径
input_files = []

# 遍历文件夹内的文件
for file in os.listdir(input_folder):
    if file.endswith(".dat"):
        file_path = os.path.join(input_folder, file)
        input_files.append(file_path)

# 对文件进行按年份排序
input_files.sort(key=lambda x: int(x[-8:-4])) # 根据你文件命名的年份位置进行调整

# 创建虚拟栅格
vrt_options = gdal.BuildVRTOptions(separate=True)  # 设置为separate以确保每个文件都在单独的波段中
vrt_file = output_file + ".vrt"
gdal.BuildVRT(vrt_file, input_files, options=vrt_options)

# 打开虚拟栅格
vrt_dataset = gdal.Open(vrt_file)

if vrt_dataset is not None:
    # 创建输出栅格
    driver = gdal.GetDriverByName("ENVI")
    output_dataset = driver.Create(output_file, vrt_dataset.RasterXSize, vrt_dataset.RasterYSize, len(input_files), gdal.GDT_Float32)

    # 将每年的数据按顺序写入对应的波段
    for i, file_path in enumerate(input_files):
        year = file_path[-8:-4]  # 从文件路径中提取年份
        band = vrt_dataset.GetRasterBand(i + 1)
        output_band = output_dataset.GetRasterBand(i + 1)
        output_band.WriteArray(band.ReadAsArray())
        output_band.SetDescription(year)  # 设置波段描述为年份
        output_band.SetNoDataValue(9999)  # 设置无效值
        output_band.FlushCache()


    # 设置地理参考信息
    output_dataset.SetGeoTransform(vrt_dataset.GetGeoTransform())
    output_dataset.SetProjection(vrt_dataset.GetProjection())

    # 关闭数据集
    vrt_dataset = None
    output_dataset = None

    os.remove(vrt_file)  # 删除虚拟栅格文件
else:
    print(f"Failed to open the virtual raster file: {vrt_file}")