ArcGIS处理Landsat9数据实战：从数据下载到可见光合成

1次阅读

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

Landsat9 是美国地质调查局（USGS）和 NASA 联合发射的最新陆地卫星，延续了 Landsat 系列长达 50 年的地球观测使命。它的 OLI- 2 传感器包含 11 个光谱波段，空间分辨率为 30 米（热红外波段为 100 米）。对于新手来说，理解这些波段特性至关重要：

波段 1 -7、9 覆盖可见光到短波红外（0.43-2.30μm）
波段 8 是全色波段（15 米分辨率）
波段 10-11 是热红外波段

可见光合成通常使用波段 4（红）、3（绿）、2（蓝）组合，这种自然色合成方案最接近人眼观察效果，适用于土地利用分类、植被监测等场景。

注册 USGS 账号（需验证邮箱）
登录后点击 ”Data Sets” 选择 ”Landsat”→”Landsat Collection 2 Level-1″
在地图界面绘制感兴趣区域或输入经纬度坐标

日期选择器设定时间范围
云量筛选建议设置 <20%（新手可放宽到 30%）
注意选择 ”Landsat 9 OLI/TIRS C2 L1″ 数据层级

优先下载整个场景（约 1GB）
压缩包包含多个波段文件（如 B2.TIF 对应蓝波段）
推荐使用 Download Manager 工具批量下载

打开 ArcGIS Pro 创建新工程
通过 Catalog 面板直接拖拽解压后的波段文件到地图窗口
或使用 ”Raster to Mosaic Dataset” 工具批量导入

注意：遇到坐标系警告时，建议统一转换为 WGS84 Web 墨卡托（EPSG:3857）

搜索打开 ”Composite Bands” 工具
按顺序添加红（B4）、绿（B3）、蓝（B2）波段
输出设置选择 ”TIFF” 格式
勾选 ”Build Pyramids” 加速后续显示

右键图层→Symbology→Stretch Type 选择 ”Percent Clip”
调整 Gamma 值（建议 1.3-1.5）增强对比度
使用 ”Enhance” 工具中的直方图均衡化

# 使用 ArcPy 进行批量投影转换
import arcpy
from arcpy.sa import *

arcpy.env.workspace = "输入文件夹路径"
raster_list = arcpy.ListRasters()

for raster in raster_list:
    out_coor_system = arcpy.SpatialReference(3857)  # Web 墨卡托
    out_raster = "输出路径 /" + raster[:-4] + "_reprojected.tif"
    arcpy.ProjectRaster_management(raster, out_raster, out_coor_system)

使用 ”Mosaic to New Raster” 工具拼接相邻场景
重叠区选择 ”Blend” 混合模式
建议保留 10% 的重叠区域

创建匀色参考点（选择亮度适中的区域）
使用 ”Colour Balance” 工具匹配直方图
或运行 ”Histogram Matching” 算法

# 完整处理流程自动化示例
import arcpy
import os

# 环境设置
arcpy.env.overwriteOutput = True
workspace = r"C:\Landsat9_Data"
output_dir = os.path.join(workspace, "Processed")
os.makedirs(output_dir, exist_ok=True)

# 1. 波段合成
def composite_bands(input_folder):
    blue = os.path.join(input_folder, "B2.TIF")
    green = os.path.join(input_folder, "B3.TIF")
    red = os.path.join(input_folder, "B4.TIF")

    out_composite = os.path.join(output_dir, "Composite.tif")
    arcpy.CompositeBands_management([red, green, blue], out_composite)
    return out_composite

# 2. 色彩增强
def enhance_image(input_raster):
    out_stretch = os.path.join(output_dir, "Enhanced.tif")
    arcpy.Stretch_management(input_raster, out_stretch, 
                           stretch_type="PERCENT_CLIP", 
                           min_percent=2, max_percent=98)
    return out_stretch

# 主流程
if __name__ == "__main__":
    scene_folder = os.path.join(workspace, "LC09_L1TP_123045_20220312")
    composite = composite_bands(scene_folder)
    final_output = enhance_image(composite)
    print(f"处理完成: {final_output}")