基于Sentinel-2的多年月度像素均值计算效率优化求助

问题背景 #

拥有6年Sentinel-2生成的NDVI、NDWI、NDMI、MNDWI、SSI、IB、IB2指数影像，需计算每个像素的多年月度均值，并导出为QGIS兼容的TIFF文件。当前代码将像素拆分为DataFrame列的方式效率极低：单张200万像素影像，仅处理10万像素列就耗时400分钟，急需优化。

原代码核心性能瓶颈 #

• 将影像像素展开为DataFrame列的操作（apply(lambda x: x.flatten())+循环创建列）是最大开销点，Pandas对大规模列操作的内存和计算成本极高
• 未利用NumPy原生的向量化运算能力，反而用循环逐像素处理，完全浪费了CPU并行计算潜力

优先方案：用NumPy/GDAL原生向量化计算（无需Dask） #

直接基于数组操作跳过DataFrame的低效转换，是最快的优化方式：

步骤1：配置参数与读取参考地理信息 #

import os
import gdal
import numpy as np

# 核心配置
indices = ['NDVI', 'NDWI', 'NDMI', 'MNDWI', 'SSI', 'IB', 'IB2']
months = ['janvier', 'fevrier', 'mars', 'avril', 'mai', 'juin', 
          'juillet', 'aout', 'sep', 'oct', 'nov', 'dec']
month_numbers = range(1,13)
root_dir = "/mnt/Data/30_Stages_Encours/2023/ZonesHumides_Sarah/sentinel2/indices"
output_dir = "/mnt/Data/30_Stages_Encours/2023/ZonesHumides_Sarah/resultats/sentinel2/cartes/qgis"

# 读取参考影像的地理参数（任意一张同规格影像即可）
ref_path = os.path.join(root_dir, "NDVI_2018_01.tif")
ds_ref = gdal.Open(ref_path)
width = ds_ref.RasterXSize
height = ds_ref.RasterYSize
gt = ds_ref.GetGeoTransform()
proj = ds_ref.GetProjection()
ds_ref = None

步骤2：批量计算月度均值并导出 #

for idx in indices:
    print(f"处理指数: {idx}")
    for month_name, month_num in zip(months, month_numbers):
        # 筛选该指数对应月份的所有年份影像
        img_paths = [os.path.join(root_dir, f) 
                     for f in os.listdir(root_dir) 
                     if idx in f and f"_{month_num:02d}.tif" in f]
        
        if not img_paths:
            print(f"无{idx}的{month_name}影像，跳过")
            continue
        
        # 批量读取影像为数组，堆叠成三维数组（年份数 × 高度 × 宽度）
        img_arrays = []
        for path in img_paths:
            ds = gdal.Open(path)
            arr = ds.GetRasterBand(1).ReadAsArray()
            img_arrays.append(arr)
            ds = None
        stack = np.stack(img_arrays, axis=0)
        
        # 计算多年月度均值（沿年份轴取平均，忽略无效值）
        mean_array = np.nanmean(stack, axis=0)
        
        # 导出为QGIS兼容的TIFF
        output_path = os.path.join(output_dir, f"{idx}_{month_name}.tif")
        driver = gdal.GetDriverByName("GTiff")
        out_ds = driver.Create(output_path, width, height, 1, gdal.GDT_Float32)
        out_ds.SetGeoTransform(gt)
        out_ds.SetProjection(proj)
        out_ds.GetRasterBand(1).WriteArray(mean_array)
        out_ds.GetRasterBand(1).SetNoDataValue(np.nan)  # 设置无效值标记
        out_ds.FlushCache()
        out_ds = None
        print(f"已生成: {output_path}")

性能提升原因 #

• 全程用NumPy向量化运算，无逐像素循环，充分利用CPU并行能力
• 跳过DataFrame列转换操作，避免不必要的内存开销和格式转换损耗
• 直接基于GDAL读写影像，中间数据处理环节最少

Dask优化方案（超大规模数据场景） #

如果影像总数据量超过内存（如100GB以上），用Dask分块并行计算避免内存溢出：

import dask.array as da
from dask_geotiff import read_geotiff

for idx in indices:
    print(f"处理指数: {idx}")
    for month_name, month_num in zip(months, month_numbers):
        img_paths = [os.path.join(root_dir, f) 
                     for f in os.listdir(root_dir) 
                     if idx in f and f"_{month_num:02d}.tif" in f]
        
        if not img_paths:
            continue
        
        # Dask分块读取所有影像（自定义分块大小适配内存）
        dask_arrays = []
        for path in img_paths:
            da_arr = read_geotiff(path, chunks=(height//4, width//4))
            dask_arrays.append(da_arr)
        stack = da.stack(dask_arrays, axis=0)
        
        # 并行计算均值（Dask自动分块处理）
        mean_da = da.nanmean(stack, axis=0)
        mean_array = mean_da.compute()  # 触发计算
        
        # 导出逻辑与之前一致
        output_path = os.path.join(output_dir, f"{idx}_{month_name}_dask.tif")
        driver = gdal.GetDriverByName("GTiff")
        out_ds = driver.Create(output_path, width, height, 1, gdal.GDT_Float32)
        out_ds.SetGeoTransform(gt)
        out_ds.SetProjection(proj)
        out_ds.GetRasterBand(1).WriteArray(mean_array)
        out_ds.GetRasterBand(1).SetNoDataValue(np.nan)
        out_ds.FlushCache()
        out_ds = None

Dask核心优势 #

• 自动分块处理数据，避免内存溢出
• 支持多线程/多进程并行，充分利用多核CPU
• 适配TB级遥感影像数据处理场景

额外优化建议 #

• 提前用gdalwarp统一所有影像的分辨率、投影、范围，避免堆叠计算出错
• 若影像存在云污染等无效值，必须用np.nanmean替代np.mean，确保均值准确性
• 若需自定义运算，可结合numba加速，但基础均值运算用NumPy/Dask已足够高效

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