批量处理geoTIFF时for循环仅导出最后一个数据集问题

脚本核心功能 #

• 遍历影像存储文件夹，识别拍摄日期相同、数量≥2张的影像组
• 将所有影像像素值插值到标准化1000m分辨率格网，修正不同影像间像素尺寸微差、像素偏移问题
• 对无重叠的同日期影像做边缘拼接
• 将处理完成的影像以GeoTIFF格式导出到指定目标文件夹

核心代码段 #

1000m分辨率标准化格网构建 #

# 创建1000m像素分辨率的标准化格网
# 计算1000m像素对应的经纬度度数
a = 6378137  # WGS84椭球长半轴（纬度方向），单位：米
b = 6356752  # WGS84椭球短半轴（经度方向），单位：米

earth_circum_lat = 2 * math.pi * a  # 地球纬度方向周长
earth_circum_lon = 2 * math.pi * b  # 地球经度方向周长

deg_lat = 360 / earth_circum_lat * 1000
deg_lon = 360 / earth_circum_lon * 1000

# 生成格网的x、y坐标序列
y = np.arange(5., -5., -deg_lat)
x = np.arange(-95., -85., deg_lon)

影像读取与日期统计 #

# 读取所有MODIS GeoTIFF文件
folder = "E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/batch"
list_of_paths = glob.glob(folder + '/*.tif', recursive=True)

# 提取文件夹内所有影像对应的拍摄日期
modis = []  # 初始化存储日期的空列表
for i in range(0, np.size(list_of_paths)):
    # 文件命名规则：cloud_effective_radius_YYYYMMDD_HHMMSS.tif
    modis.append(list_of_paths[i].split('ius_')[1][0:8])

# 去重得到文件夹内所有不同的拍摄日期
modis = np.unique(modis)
print(modis)
print('\n该文件夹内共包含{}个不同日期的数据\n'.format(np.size(modis)))

重格网、拼接与导出逻辑 #

# GeoTIFF插值、拼接批量处理主循环
for i in range(0, np.size(modis)):
    # 提取当前日期对应的所有影像路径
    list_files_date = glob.glob(os.path.join(folder, 'cloud_effective_radius_{0}*.tif'.format(modis[i])))

    # 当日影像数量大于1时执行拼接
    if len(list_files_date) > 1:
        ds = rioxarray.open_rasterio(list_files_date[0], engine='rasterio')
        ds1 = rioxarray.open_rasterio(list_files_date[1], engine='rasterio')

        # 重格网：将所有影像插值到统一坐标格网
        ds_interp = ds.interp(y=y, x=x, method="nearest")
        ds1_interp = ds1.interp(y=y, x=x, method="nearest")

        # 逐变量合并数据集
        ds_merged = ds_interp.combine_first(ds1_interp)
        # 合并规则：两数据集同位置均为非空值时保留第一个的值；仅一个为空时取非空值

        if np.size(list_files_date) > 2:
            for j in range(2, np.size(list_files_date)):
                ds = rioxarray.open_rasterio(list_files_date[j], engine='rasterio')
                ds_interp = ds.interp(y=y, x=x, method="nearest")
                ds_merged = ds_merged.combine_first(ds_interp)

    # 当日仅1张影像时仅执行重格网
    else:
        list_files_date = (glob.glob(os.path.join(folder, 'cloud_effective_radius_{0}*.tif'.format(modis[i]))))
        ds = rioxarray.open_rasterio(list_files_date[0], engine='rasterio')
        ds_merged = ds.interp(y=y, x=x, method="nearest")

    # 导出为GeoTIFF文件
    img_number = 1
    ds_merged.rio.to_raster("E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/2021_06_merged/ds_merged" +
                            str(img_number) + ".tif",
                            driver="GTiff")
    img_number += 1

问题现象 #

运行脚本后目标文件夹内仅生成1个合并后的GeoTIFF文件，始终对应最后一个日期的处理结果，推测前期处理结果被覆盖，但无法定位具体原因：原逻辑中img_number += 1设计目的是避免文件重名覆盖，无法判断问题是源于导出代码逻辑错误，还是for循环内ds_merged变量被重复覆盖导致。

问题根因 #

核心问题是计数器变量初始化位置错误

• img_number = 1被写在了日期遍历的for循环内部，每次循环启动时该变量都会被重置为1，导致每次循环导出的文件名固定为ds_merged1.tif，前一次循环生成的文件会被后一次的写入操作直接覆盖，最终仅保留最后一个日期的处理结果。
• for循环内ds_merged变量重复赋值属于正常逻辑：每次循环处理单个日期的数据，处理完成后立即执行导出，不存在变量覆盖导致前期数据丢失的问题。

修复方案 #

方案1：调整计数器初始化位置 #

将img_number = 1移动到for循环外部，保证每次循环导出的文件名序号递增，避免重名覆盖：

# 计数器初始化移到循环外部
img_number = 1
for i in range(0, np.size(modis)):
    # 原有插值、拼接逻辑保持不变
    
    # 导出文件
    ds_merged.rio.to_raster("E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/2021_06_merged/ds_merged" +
                            str(img_number) + ".tif",
                            driver="GTiff")
    img_number += 1

方案2：用日期作为文件名（推荐） #

序号命名的文件无法直接对应拍摄日期，检索效率低，可直接用提取到的日期字符串作为文件名，从根源上避免重名，也方便后续数据管理：

for i in range(0, np.size(modis)):
    current_date = modis[i]
    # 原有插值、拼接逻辑保持不变
    
    # 以日期命名导出文件
    output_path = f"E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/2021_06_merged/ds_merged_{current_date}.tif"
    ds_merged.rio.to_raster(output_path, driver="GTiff")

内容的提问来源于stack exchange，提问作者captainjasper