Python卫星影像处理求助：每2个月重采样并计算波段均值

1. 解析文件名并构建元数据DataFrame #

先批量提取所有影像的关键信息（tile、日期、波段、文件路径），用Pandas构建DataFrame，方便后续分组和时间处理：

import os
import pandas as pd
from datetime import datetime

# 替换为你的影像存储路径
img_dir = "satellite_images"
file_list = [f for f in os.listdir(img_dir) if f.endswith(".tif")]

metadata = []
for fname in file_list:
    parts = fname.split("_")
    tile = parts[0]
    # 提取日期部分（去掉时间戳）并转为datetime格式
    date_str = parts[1].split("T")[0]
    date = datetime.strptime(date_str, "%Y%m%d")
    band = parts[2]
    metadata.append({
        "tile": tile,
        "date": date,
        "band": band,
        "file_path": os.path.join(img_dir, fname)
    })

# 将元数据转为DataFrame并设置日期为索引
df = pd.DataFrame(metadata)
df.set_index("date", inplace=True)

2. 按Tile+波段分组，执行2个月时间重采样与均值计算 #

通过groupby按tile和波段分组，对每组的时间序列做2个月重采样，再计算窗口内的影像均值：

import rasterio
import numpy as np

def calculate_band_mean(group):
    # 读取组内所有影像数据，保留原影像的地理信息
    arrays = []
    profile = None
    for _, row in group.iterrows():
        with rasterio.open(row["file_path"]) as src:
            if not profile:
                profile = src.profile.copy()
            # 读取单波段数据，同时将nodata转为NaN
            arr = src.read(1)
            arr[arr == src.nodata] = np.nan
            arrays.append(arr)
    
    # 计算时间维度的均值（忽略NaN）
    mean_arr = np.nanmean(np.stack(arrays, axis=0), axis=0)
    
    # 更新输出影像的配置（比如将数据类型转为float32存储均值）
    profile.update(dtype=rasterio.float32, count=1)
    
    # 生成输出文件名，示例：T29PQN_B01_202205_mean.tif
    tile = group["tile"].iloc[0]
    band = group["band"].iloc[0]
    window_start = group.index[0].strftime("%Y%m")
    output_path = os.path.join("output_mean_images", f"{tile}_{band}_{window_start}_mean.tif")
    
    # 写入均值影像
    with rasterio.open(output_path, "w", **profile) as dst:
        dst.write(mean_arr, 1)
    
    return mean_arr

# 按tile和band分组，每2个月重采样（"2MS"表示以每2个月的第一天为窗口起始，可按需改为"2M"）
df.groupby(["tile", "band"]).resample("2MS").apply(calculate_band_mean)

3. 关键优化与注意事项 #

• 频率参数调整：Pandas重采样频率"2MS"是每2个月的月初，"2M"是月末，可根据你的时间序列需求切换。
• 内存优化：如果影像分辨率高、数量多，可改用rasterio的分块读取（src.block_windows()），避免一次性加载所有数据占满内存。
• 并行加速：1836景数据量较大，可借助multiprocessing或dask对分组任务并行处理，大幅缩短运行时间。
• 掩膜校验：确认原影像的nodata值设置正确，计算均值时用np.nanmean确保掩膜区域不参与计算。

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