共计 3575 个字符，预计需要花费 9 分钟才能阅读完成。
背景介绍
Landsat 系列卫星作为历史最悠久的对地观测计划之一，提供了连续 50 年的地球表面数据。其多光谱数据具有 30 米空间分辨率（全色波段 15 米）和 16 天重访周期，广泛应用于土地利用分类、植被监测、环境变化分析等领域。对于 GIS 新手而言，掌握 Landsat 数据处理是开展遥感分析的基础技能。
技术选型对比
- ArcGIS 优势：
- 无缝集成空间分析与制图功能
- 可视化操作界面降低学习门槛
- ArcPy 模块支持 Python 脚本自动化
-
企业级数据管理能力
-
ENVI 特点：
- 专业的辐射校正工具链
- 更丰富的光谱分析功能
-
适合科研级深度处理
-
QGIS 优势：
- 开源免费
- 丰富的插件生态系统
- 轻量化处理中小规模数据
核心实现步骤
数据获取
- 访问 USGS EarthExplorer（https://earthexplorer.usgs.gov/）
- 注册账号后选择 Landsat Collection 2 Level- 1 数据
- 按以下条件筛选：
- 云量覆盖率 <20%
- 包含所有标准波段（B1-B11）
- 选择 UTM 投影的 GeoTIFF 格式
预处理流程
辐射定标
将 DN 值转换为大气顶层反射率（TOA）：
import arcpy
from arcpy.sa import *
# 设置工作空间
arcpy.env.workspace = r"D:\Landsat\LC08_L1TP_123032_20220101"
# 辐射定标参数（以 Landsat 8 为例）band_names = ["B2", "B3", "B4"] # 蓝、绿、红波段
refl_mult = 0.0000275 # 反射率乘数
refl_add = -0.2 # 反射率加数
# 批量处理波段
for band in band_names:
in_raster = band + ".TIF"
out_raster = band + "_TOA"
arcpy.gp.RasterCalculator_sa(f"{refl_mult} * Float({in_raster}) + {refl_add}", out_raster)大气校正（快速黑暗像元法）
# 假设已获取暗像元值（通过 histogram 分析）dark_value = 0.01
for band in band_names:
toa_raster = band + "_TOA"
corrected = (Raster(toa_raster) - dark_value) / (1 - dark_value)
corrected.save(band + "_CORRECTED")波段组合
真彩色合成（RGB=4,3,2）
# 使用 CompositeBands 工具
input_rasters = ["B4_CORRECTED", "B3_CORRECTED", "B2_CORRECTED"]
arcpy.CompositeBands_management(input_rasters, "TrueColor.tif")假彩色合成（植被增强 RGB=5,4,3）
input_rasters = ["B5_CORRECTED", "B4_CORRECTED", "B3_CORRECTED"]
arcpy.CompositeBands_management(input_rasters, "FalseColor.tif")完整代码示例
import arcpy
import logging
from datetime import datetime
# 配置日志
logging.basicConfig(filename='landsat_processing.log',
level=logging.INFO,
format='%(asctime)s - %(levelname)s - %(message)s')
def process_landsat(input_folder, output_folder):
try:
arcpy.env.workspace = input_folder
arcpy.env.overwriteOutput = True
# 1. 辐射定标
logging.info("Starting radiometric calibration")
bands = {"B2": "Blue", "B3": "Green", "B4": "Red", "B5": "NIR"}
for band in bands.keys():
in_raster = f"{band}.TIF"
if arcpy.Exists(in_raster):
out_toa = f"{output_folder}\\{band}_TOA.tif"
arcpy.gp.RasterCalculator_sa("0.0000275 * Float(\"" + in_raster + "\") - 0.2",
out_toa)
# 2. 大气校正
logging.info("Performing atmospheric correction")
dark_value = 0.01
for band in bands.keys():
toa_raster = f"{output_folder}\\{band}_TOA.tif"
if arcpy.Exists(toa_raster):
corrected = (Raster(toa_raster) - dark_value) / (1 - dark_value)
corrected.save(f"{output_folder}\\{band}_CORRECTED.tif")
# 3. 波段合成
logging.info("Creating band composites")
true_color = [f"{output_folder}\\B4_CORRECTED.tif",
f"{output_folder}\\B3_CORRECTED.tif",
f"{output_folder}\\B2_CORRECTED.tif"]
if all(arcpy.Exists(r) for r in true_color):
arcpy.CompositeBands_management(
true_color,
f"{output_folder}\\TrueColor_Composite.tif")
# 添加金字塔和统计信息
logging.info("Building pyramids")
arcpy.BuildPyramids_management(output_folder)
logging.info("Processing completed successfully")
return True
except Exception as e:
logging.error(f"Processing failed: {str(e)}", exc_info=True)
return False
# 使用示例
if __name__ == "__main__":
input_path = r"D:\Landsat\LC08_L1TP_123032_20220101"
output_path = r"D:\Processed\Output"
process_landsat(input_path, output_path)性能优化技巧
-
分块处理：
arcpy.env.tileSize = "256 256" arcpy.env.compression = "LZW" -
内存管理：
- 设置临时工作空间到 SSD
- 处理前使用
arcpy.Delete_management()清理临时数据 -
分幅处理大区域数据
-
并行计算：
arcpy.env.parallelProcessingFactor = "75%"
常见问题解决
-
投影错误：
# 强制统一投影 arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(32651) # WGS84 UTM 51N -
波段顺序混乱：
- 使用
arcpy.Describe()检查波段属性 -
重建波段列表时按波长排序
-
色偏问题：
- 检查辐射定标系数是否正确
- 应用直方图均衡化增强：
arcpy.Stretch_management("Composite.tif", "Enhanced.tif", "PercentClip", "0 1", "0 255", "NONE")
实际应用案例
- 农业监测：
- 结合 NDVI 识别作物健康状况
-
时序分析监测生长周期
-
城市扩张：
- 多期影像对比检测建设用地变化
-
使用假彩色合成突出建成区
-
灾害评估：
- 洪涝灾害前后影像对比
- 使用 NDWI 提取水体范围
延伸思考
- 将处理流程封装为 ArcGIS Toolbox 工具
- 结合深度学习进行地物分类
- 搭建自动化处理服务（使用 ArcGIS Server 或 ArcGIS Pro Notebooks）
结语
通过本文介绍的完整流程，GIS 开发者可以系统掌握 Landsat 数据在 ArcGIS 平台的处理方法。建议新手先从单景影像入手，逐步扩展到时序分析和批量处理。实际应用中需注意不同 Landsat 传感器（如 OLI 与 TM）的参数差异，持续积累处理经验才能高效应对各种业务场景需求。
正文完
