使用YATSM时Landsat影像堆叠报错求助，求Rasterio/Rio替代方案

1. 同类YATSM预处理问题的常见情况 #

• 关于batch_landsat模块缺失：这个模块是早期YATSM版本的配套工具，后来的版本要么整合进了核心代码，要么直接移除了，很多用户在复用旧脚本时都会碰到这个问题。如果非要用原生YATSM工具链，你可以尝试安装对应旧版本的YATSM依赖；但更高效的方式是直接换用通用栅格处理工具（比如下面的Rasterio方案），避免版本兼容坑。
• 关于TypeError: 'NoneType' object is not subscriptable：这个错误几乎都是因为landsat_stack.py没正确读取到影像的范围（extent），大概率是你的Landsat影像元数据（比如MTL文件）丢失、影像投影不一致，或者脚本没正确关联MTL文件来获取范围参数。不少用户反馈过类似问题，排查方向优先检查MTL文件是否存在、输入路径是否正确、所有影像是否处于同一投影。

2. 基于Rasterio的替代解决方案 #

下面是一个完全自定义的Python脚本，用Rasterio实现你需要的8波段堆叠，严格匹配你的格式要求：

前提准备 #

先确保安装好依赖：

pip install rasterio numpy glob2

堆叠脚本（适配Landsat 8/9，可按需调整） #

这个脚本会遍历指定目录下的Landsat景，按要求处理每个波段后堆叠成单文件，你可以根据自己的文件命名规则微调波段匹配逻辑：

import os
import numpy as np
import rasterio
import glob2
from rasterio.warp import calculate_default_transform, reproject, Resampling

# -------------------------- 配置参数 --------------------------
input_root = "/path/to/your/landsat/scenes"  # 存放所有Landsat景的根目录
output_root = "/path/to/your/output/stacked" # 堆叠影像输出目录
target_crs = "EPSG:32650"                    # 统一目标投影（替换为你的UTM带EPSG码）
# -------------------------------------------------------------

os.makedirs(output_root, exist_ok=True)

# 遍历每个Landsat景目录
for scene_name in os.listdir(input_root):
    scene_path = os.path.join(input_root, scene_name)
    if not os.path.isdir(scene_path):
        continue

    # 匹配各个波段文件（根据你的文件名格式调整通配符）
    band_mapping = {
        "sr_b1": glob2.glob(os.path.join(scene_path, "*SR_B1.TIF"))[0],
        "sr_b2": glob2.glob(os.path.join(scene_path, "*SR_B2.TIF"))[0],
        "sr_b3": glob2.glob(os.path.join(scene_path, "*SR_B3.TIF"))[0],
        "sr_b4": glob2.glob(os.path.join(scene_path, "*SR_B4.TIF"))[0],
        "sr_b5": glob2.glob(os.path.join(scene_path, "*SR_B5.TIF"))[0],
        "sr_b7": glob2.glob(os.path.join(scene_path, "*SR_B7.TIF"))[0],
        "thermal_b6": glob2.glob(os.path.join(scene_path, "*ST_B6.TIF"))[0],
        "fmask": glob2.glob(os.path.join(scene_path, "*Fmask.TIF"))[0]
    }

    processed_bands = []
    output_profile = None

    # 处理SR波段：乘10000，转uint16
    for sr_key in ["sr_b1", "sr_b2", "sr_b3", "sr_b4", "sr_b5", "sr_b7"]:
        with rasterio.open(band_mapping[sr_key]) as src:
            if output_profile is None:
                # 以第一个SR波段为基础构建输出profile
                output_profile = src.profile.copy()
                output_profile.update(count=8, dtype="uint16")
            # 读取并转换数据
            data = src.read(1)
            data = (data * 10000).astype(np.uint16)
            # 统一投影（如果当前波段投影不匹配）
            if src.crs != target_crs:
                transform, width, height = calculate_default_transform(
                    src.crs, target_crs, src.width, src.height, *src.bounds
                )
                reproj_data = np.zeros((height, width), dtype=np.uint16)
                reproject(
                    source=data,
                    destination=reproj_data,
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=target_crs,
                    resampling=Resampling.bilinear
                )
                processed_bands.append(reproj_data)
            else:
                processed_bands.append(data)

    # 处理热红外波段：乘100，转uint16
    with rasterio.open(band_mapping["thermal_b6"]) as src:
        data = src.read(1)
        data = (data * 100).astype(np.uint16)
        if src.crs != target_crs:
            transform, width, height = calculate_default_transform(
                src.crs, target_crs, src.width, src.height, *src.bounds
            )
            reproj_data = np.zeros((height, width), dtype=np.uint16)
            reproject(
                source=data,
                destination=reproj_data,
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=target_crs,
                resampling=Resampling.bilinear
            )
            processed_bands.append(reproj_data)
        else:
            processed_bands.append(data)

    # 处理Fmask：直接读取，保持uint8类型
    with rasterio.open(band_mapping["fmask"]) as src:
        data = src.read(1).astype(np.uint8)
        if src.crs != target_crs:
            transform, width, height = calculate_default_transform(
                src.crs, target_crs, src.width, src.height, *src.bounds
            )
            reproj_data = np.zeros((height, width), dtype=np.uint8)
            reproject(
                source=data,
                destination=reproj_data,
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=target_crs,
                resampling=Resampling.nearest  # 分类数据用最近邻重采样
            )
            processed_bands.append(reproj_data)
        else:
            processed_bands.append(data)

    # 写入堆叠后的影像
    output_path = os.path.join(output_root, f"{scene_name}_stacked.tif")
    with rasterio.open(output_path, "w", **output_profile) as dst:
        for idx, band in enumerate(processed_bands, 1):
            dst.write(band, idx)
    
    print(f"✅ 完成景 {scene_name} 的堆叠")

脚本调整要点 #

• 文件匹配：如果是Landsat 5/7，热红外波段是B6而非ST_B6，Fmask文件名可能是小写的*fmask.tif，需要对应修改band_mapping里的通配符。
• 投影设置：把target_crs替换为你需要的UTM投影EPSG码，确保所有影像统一到同一投影，彻底避免范围不匹配问题。
• 数据类型：SR波段乘10000后用uint16足够存储（0-10000对应0-1的反射率），热红外乘100后用uint16能覆盖常规温度范围，Fmask用uint8正好匹配0-255的分类值。

批量调用的bash脚本 #

如果习惯用bash触发批量处理，可以写一个简单的run_stack.sh：

#!/bin/bash
INPUT_DIR="/path/to/your/landsat/scenes"
OUTPUT_DIR="/path/to/your/output/stacked"
TARGET_CRS="EPSG:32650"

python3 stack_landsat.py --input-dir $INPUT_DIR --output-dir $OUTPUT_DIR --target-crs $TARGET_CRS

（需要把Python脚本改成支持argparse命令行参数的版本，这样更灵活）

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