使用 SNAP-Python（snappy）处理 Sentinel-1：安装与批处理

之前在做积雪深度反演项目的时候，需要批量处理 Sentinel-1 数据。由于研究区范围较大，手工操作效率较低，于是选择使用 SNAP-Python（snappy）构建自动化预处理流程。

SNAP-Python 是 ESA SNAP 软件的 Python 接口，能够直接调用 SNAP 算子处理 Sentinel 数据，同时可以与 Python 生态协同（如 NumPy、Matplotlib 等），非常适合搭建可复现、可扩展的遥感处理脚本。

导读

本文包含三部分内容：

1. 安装入口与关键注意点
2. 一套可直接参考的 Sentinel-1 预处理代码
3. 常见报错排查与参数速查表

处理流程总览

flowchart TD
        A[读取 Sentinel-1 产品] --> B[应用精确轨道文件]
        B --> C[热噪声去除]
        C --> D[辐射定标]
        D --> E[斑点滤波]
        E --> F[地形校正]
        F --> G[按研究区裁剪]
        G --> H[导出 GeoTIFF]

安装

首先安装 SNAP-Python，可参考官方文档。

关键步骤可参考下图：

使用 SNAP-Python 处理 Sentinel-1 数据

安装完成后，就可以开始处理 Sentinel-1 数据。下面给出一个完整示例，演示如何读取 Sentinel-1 产品并执行标准预处理流程。

▶

使用SNAP-Python处理Sentinel-1数据

# 示例代码：使用SNAP-Python处理Sentinel-1数据
# 读取解压后的Sentinel-1 GRD产品(EW和IW模式)

import datetime
import time
from esa_snappy import ProductIO  # SNAP的产品输入输出模块
from esa_snappy import HashMap   # SNAP的参数映射容器
import os, gc                    # 操作系统接口和垃圾回收
from esa_snappy import GPF       # SNAP的图像处理框架(Graph Processing Framework)

def do_apply_orbit_file(source):
    """
    应用精确轨道文件
    精确轨道文件用于提高卫星位置精度，改善几何校正质量
    
    参数:
        source: 输入的SAR产品
    返回:
        output: 应用轨道文件后的产品
    """
    print('\tApply orbit file...')
    parameters = HashMap()
    parameters.put('Apply-Orbit-File', True)  # 启用轨道文件应用
    output = GPF.createProduct('Apply-Orbit-File', parameters, source)
    return output

def do_thermal_noise_removal(source):
    """
    热噪声去除
    去除Sentinel-1数据中的热噪声，特别是在图像边缘和暗区域
    
    参数:
        source: 输入的SAR产品
    返回:
        output: 去除热噪声后的产品
    """
    print('\tThermal noise removal...')
    parameters = HashMap()
    parameters.put('removeThermalNoise', True)  # 启用热噪声去除
    output = GPF.createProduct('ThermalNoiseRemoval', parameters, source)
    return output

def do_calibration(source, polarization, pols):
    """
    辐射定标
    将数字值转换为物理量(如后向散射系数sigma0)
    
    参数:
        source: 输入的SAR产品
        polarization: 极化模式(DH/DV/SH/SV等)
        pols: 极化字符串(如'VH,VV')
    返回:
        output: 定标后的产品
    """
    print('\tCalibration...')
    parameters = HashMap()
    parameters.put('outputSigmaBand', True)  # 输出sigma0波段
    
    # 根据不同极化模式选择相应的源波段
    if polarization == 'DH':                 # 双极化HH+HV
        parameters.put('sourceBands', 'Intensity_HH,Intensity_HV')
    elif polarization == 'DV':              # 双极化VH+VV
        parameters.put('sourceBands', 'Intensity_VH,Intensity_VV')
    elif polarization == 'SH' or polarization == 'HH':  # 单极化HH
        parameters.put('sourceBands', 'Intensity_HH')
    elif polarization == 'SV':              # 单极化VV
        parameters.put('sourceBands', 'Intensity_VV')
    else:
        print("different polarization!")
    
    parameters.put('selectedPolarisations', pols)      # 选择的极化
    parameters.put('outputImageScaleInDb', False)      # 输出线性尺度(不转换为dB)
    output = GPF.createProduct("Calibration", parameters, source)
    return output

def do_speckle_filtering(source):
    """
    斑点噪声滤波
    使用Lee滤波器去除SAR图像中的斑点噪声，提高图像质量
    
    参数:
        source: 输入的SAR产品
    返回:
        output: 滤波后的产品
    """
    print('\tSpeckle filtering...')
    parameters = HashMap()
    parameters.put('filter', 'Lee')          # 使用Lee滤波器
    parameters.put('SizeX', 5)         # X方向滤波窗口大小为5像素
    parameters.put('SizeY', 5)         # Y方向滤波窗口大小为5像素
    output = GPF.createProduct('Speckle-Filter', parameters, source)
    return output

def do_terrain_correction(source, proj, downsample):
    """
    地形校正
    使用数字高程模型(DEM)进行地形校正，消除地形起伏造成的几何畸变
    
    参数:
        source: 输入的SAR产品
        proj: 地图投影参数字符串
        downsample: 是否需要降采样(1-需要降采样到40m，0-不需要)
    返回:
        output: 地形校正后的产品
    """
    print('\tTerrain correction...')
    parameters = HashMap()
    parameters.put('demName', 'GETASSE30')                    # 使用GETASSE30 DEM
    parameters.put('imgResamplingMethod', 'BILINEAR_INTERPOLATION')  # 双线性插值重采样
    parameters.put('mapProjection', proj)                     # 设置地图投影(如UTM)，注释后为默认WGS84
    parameters.put('saveProjectedLocalIncidenceAngle', True)  # 保存投影后的局地入射角
    parameters.put('saveSelectedSourceBand', True)            # 保存选中的源波段
    
    # 如果需要降采样，设置像素间距为40米
    while downsample == 1:                                    
        parameters.put('pixelSpacingInMeter', 40.0)           # 设置像素间距为40米
        break
    output = GPF.createProduct('Terrain-Correction', parameters, source)
    return output

def do_subset(source, wkt):
    """
    子集裁剪
    根据指定的地理区域裁剪图像
    
    参数:
        source: 输入的SAR产品
        wkt: Well-Known Text格式的地理区域定义
    返回:
        output: 裁剪后的产品
    """
    print('\tSubsetting...')
    parameters = HashMap()
    parameters.put('geoRegion', wkt)          # 设置地理区域
    output = GPF.createProduct('Subset', parameters, source)
    return output

def main():
    """
    主函数 - 批量处理Sentinel-1数据的完整预处理流程
    
    处理步骤:
    1. 应用精确轨道文件
    2. 热噪声去除
    3. 辐射定标
    4. 斑点滤波
    5. 地形校正
    6. 子集裁剪
    7. 输出为GeoTIFF格式
    """
    ## 设置输入和输出路径
    path = r'data'              # Sentinel-1数据的输入路径
    outpath = r's1_preprocessed'     # 预处理结果的输出路径
    
    # 如果输出目录不存在，则创建
    if not os.path.exists(outpath):
        os.makedirs(outpath)
    
    ## 定义研究区域的WKT多边形(可通过SNAP软件绘制多边形获取)
    wkt = 'POLYGON((112.7083 28.937,113.4886 28.937,113.4886 29.5445,112.7083 29.5445,112.7083 28.937))'  # 示例区域(湖南省岳阳市)
      
    ## 投影设置使用参考arcgis中的投影参数
    proj = 'EPSG:32649'  # UTM第49区投影参数(适用于湖南岳阳地区)
    # proj = 'EPSG:4326'  # WGS84地理坐标系(如果不需要投影转换) 

    # 遍历输入路径中的所有.SAFE文件夹
    for folder in os.listdir(path):
        gc.enable()                        # 启用垃圾回收
        gc.collect()                       # 执行垃圾回收，释放内存
        
        # 读取Sentinel-1产品
        sentinel_1 = ProductIO.readProduct(path + "\\" + folder + "\\manifest.safe")
        print(sentinel_1)

        # 记录处理开始时间
        loopstarttime=str(datetime.datetime.now())
        print('Start time:', loopstarttime)
        start_time = time.time()

        ## 从文件名中提取模式、产品类型和极化信息
        # Sentinel-1文件名格式: S1A_MODE_TYPE_POLARIZATION_...
        modestamp = folder.split("_")[1]    # 获取模式(IW/EW)
        productstamp = folder.split("_")[2] # 获取产品类型(GRDH/GRDM)
        polstamp = folder.split("_")[3]     # 获取极化信息

        # 解析极化模式
        polarization = polstamp[2:4]        # 提取极化标识
        if polarization == 'DV':            # 双极化VH+VV
            pols = 'VH,VV'
        elif polarization == 'DH':          # 双极化HH+HV
            pols = 'HH,HV'
        elif polarization == 'SH' or polarization == 'HH':  # 单极化HH
            pols = 'HH'
        elif polarization == 'SV':          # 单极化VV
            pols = 'VV'
        else:
            print("Polarization error!")

        ## 开始预处理流程
        print("Starting preprocessing pipeline...")
        
        # 步骤1: 应用精确轨道文件
        applyorbit = do_apply_orbit_file(sentinel_1)
        
        # 步骤2: 热噪声去除
        thermaremoved = do_thermal_noise_removal(applyorbit)
        
        # 步骤3: 辐射定标
        calibrated = do_calibration(thermaremoved, polarization, pols)
        
        # 步骤4: 斑点滤波
        down_filtered = do_speckle_filtering(calibrated)
        
        # 释放中间产品的内存
        del applyorbit
        del thermaremoved
        del calibrated
        
        ## 根据模式和产品类型进行不同的处理
        # IW GRDH和EW GRDH产品: 降采样到40m分辨率
        if (modestamp == 'IW' and productstamp == 'GRDH') or (modestamp == 'EW' and productstamp == 'GRDH'):
            # 步骤5: 地形校正(带降采样)
            down_tercorrected = do_terrain_correction(down_filtered, proj, 1)
            # 步骤6: 子集裁剪
            down_subset = do_subset(down_tercorrected, wkt)
            
            # 释放内存
            del down_filtered
            del down_tercorrected
            
        # EW GRDM产品: 保持原始分辨率(约40m)
        elif modestamp == 'EW' and productstamp == 'GRDM':
            # 步骤5: 地形校正(无降采样)
            tercorrected = do_terrain_correction(down_filtered, proj, 0)
            # 步骤6: 子集裁剪
            subset = do_subset(tercorrected, wkt)
            
            # 释放内存
            del down_filtered
            del tercorrected
        else:
            print("Different spatial resolution is found.")

        ## 确定输出文件并写入结果
        down = 1
        try: 
            down_subset                      # 检查是否存在降采样结果
        except NameError:
            down = None                      # 如果不存在，设置为None
            
        if down is None:
            # 输出原始分辨率结果
            print("Writing...")
            ProductIO.writeProduct(subset, outpath + '\\' + folder[:-5], 'GeoTIFF')
            del subset
        elif down == 1:
            # 输出降采样结果(文件名添加'_40'标识)
            print("Writing undersampled image...")
            ProductIO.writeProduct(down_subset, outpath + '\\' + folder[:-5] + '_40', 'GeoTIFF')
            del down_subset
        else:
            print("Error.")

        # 清理资源
        print('Done.')
        sentinel_1.dispose()                 # 释放产品资源
        sentinel_1.closeIO()                 # 关闭IO流
        print("--- %s seconds ---" % (time.time() - start_time))  # 显示处理耗时

if __name__== "__main__":
    main()

常见问题排查

1. 提示找不到 esa_snappy
   • 检查 SNAP 与 Python 版本是否匹配。
   • 按官方文档重新执行接口配置步骤。
2. 轨道文件步骤失败
   • 通常与网络或缓存有关，建议检查网络后重试。
   • 可先在 SNAP 图形界面验证同一数据是否可正常处理。
3. 内存不足或处理中断
   • 减小批处理规模，先分批跑。
   • 适当保留或加强中间结果释放逻辑，避免长时间占用内存。
4. 输出投影或分辨率不符合预期
   • 重点检查 mapProjection 与 pixelSpacingInMeter 参数。
   • 检查 IW/EW、GRDH/GRDM 分支逻辑是否命中正确。

参数速查表

参数	作用	示例值	备注
mapProjection	输出投影	EPSG:32649	可按研究区替换
pixelSpacingInMeter	输出分辨率（米）	40.0	仅在需要降采样时设置
selectedPolarisations	极化通道	VH,VV / HH,HV	与产品极化一致
demName	地形校正使用 DEM	GETASSE30	可按精度需求调整
geoRegion	裁剪区域（WKT）	POLYGON((…))	建议先在软件中验证范围

使用建议

1. 先用少量影像验证参数，再批量处理全量数据。
2. 优先确认投影、极化与裁剪区域配置，避免返工。
3. 处理大量影像时，建议记录日志与耗时，便于后续优化。

小结

这套流程适合 Sentinel-1 常规预处理，尤其适用于大范围、多景影像的自动化场景。你可以在此基础上继续扩展，例如增加多时相统计、质量控制日志或并行处理，以进一步提升处理效率与稳定性。