基于ArcGIS的Landsat数据处理实战：从多光谱到合成可见光

1次阅读

没有评论

共计 2619 个字符，预计需要花费 7 分钟才能阅读完成。

Landsat 卫星数据是遥感领域最常用的数据源之一，但直接将原始波段合成可见光图像时，经常会遇到色彩失真的问题。这主要是因为：

原始 DN 值（Digital Number）未经辐射定标，无法真实反映地表反射率
不同波段间的动态范围差异导致合成后整体色调偏暗或偏亮
大气散射等因素未校正，图像呈现雾蒙蒙的效果

传统 ENVI/Pix4D 等专业软件虽然能处理这些问题，但在批量处理大量数据时存在明显瓶颈：

操作流程繁琐，需要多次手动设置参数
处理大规模影像时内存占用高，容易崩溃
缺乏有效的自动化机制，重复工作量大

ArcGIS Pro 的 Composite Bands 工具是处理多光谱数据的基础。关键参数设置如下：

输入波段顺序必须对应 RGB 通道。对于 Landsat 8/9：
红色通道：Band 4（0.64–0.67 μm）
绿色通道：Band 3（0.53–0.59 μm）
蓝色通道：Band 2（0.45–0.51 μm）
输出像元类型建议选择 ”32_BIT_FLOAT” 保留精度

使用 Raster Function 模板实现 DN 值转反射率的公式：

# Landsat 8/ 9 辐射定标公式（反射率）reflectance = (DN * 0.00002) - 0.1

完整 Python 脚本示例（包含异常处理和日志）：

import arcpy
import os
import logging

# 配置日志
logging.basicConfig(filename='landsat_processing.log', level=logging.INFO)

try:
    # 输入输出路径
    input_folder = r"C:\Data\Landsat"
    output_folder = r"C:\Output"

    # 获取所有 TIFF 文件
    tif_files = [f for f in os.listdir(input_folder) if f.endswith('.TIF')]

    # 按波段类型分类
    bands = {}
    for f in tif_files:
        if '_B2.' in f: bands['blue'] = os.path.join(input_folder, f)
        elif '_B3.' in f: bands['green'] = os.path.join(input_folder, f)
        elif '_B4.' in f: bands['red'] = os.path.join(input_folder, f)

    # 检查波段完整性
    if len(bands) != 3:
        raise ValueError("缺失必要波段文件")

    # 波段合成
    logging.info("开始波段合成处理")
    composite_output = os.path.join(output_folder, "composite.tif")
    arcpy.management.CompositeBands([bands['red'], bands['green'], bands['blue']],
        composite_output
    )

    # 辐射定标
    def apply_reflectance(in_raster):
        expr = "(Float(Raster(\"{}\")) * 0.00002) - 0.1".format(in_raster)
        out_raster = arcpy.sa.RasterCalculator([in_raster], ["x"], expr, "32_BIT_FLOAT"
        )
        return out_raster

    calibrated_output = os.path.join(output_folder, "calibrated.tif")
    apply_reflectance(composite_output).save(calibrated_output)

    logging.info("处理完成")

except Exception as e:
    logging.error(f"处理失败: {str(e)}")
    raise

当处理大区域影像时，不同投影转换可能导致像元错位。解决方案：

在处理前统一所有波段的坐标系
使用 ”NEAREST” 重采样方法保持像元值不变
检查像元对齐状态：

arcpy.management.AlignRasters(["band2.tif", "band3.tif", "band4.tif"],
    "aligned_band2.tif",
    resample_type="NEAREST"
)

Collection 1 和 Collection 2 数据的主要区别：

Collection 2 使用改进的大气校正算法
元数据文件结构不同（MTL.txt 格式变化）
建议统一使用 Collection 2 数据，处理前检查元数据版本

处理大范围影像时的优化方法：

启用金字塔构建加速显示
使用分块处理（Tile Size 建议 512×512）
关闭不必要的背景处理：

arcpy.env.compression = "LZW"
arcpy.env.pyramid = "PYRAMIDS -1 BILINEAR DEFAULT 75 NO_SKIP"
arcpy.env.tileSize = "512 512"

使用直方图匹配确保多景影像色彩一致：

reference = r"C:\BaseImage.tif"
target = r"C:\NewImage.tif"
matched_output = arcpy.ia.HistogramMatch(target, reference)
matched_output.save(r"C:\Matched.tif")

两种主流压缩算法的性能比较：

压缩算法	压缩率	写入速度	CPU 占用
LZW	中等	较快	低
DEFLATE	高	慢	高

建议选择：
– 存储空间紧张时用 DEFLATE
– 需要快速读写时用 LZW

功能	ArcGIS Pro 方案	QGIS 方案
波段合成	Composite Bands 工具	波段计算器
辐射定标	Raster Function 模板	GRASS r.mapcalc
自动化处理	arcpy 脚本	PyQGIS 脚本
内存管理	分块处理参数完善	依赖系统设置
色彩校正	Histogram Match 函数	直方图拉伸插件