通过 Snappy Python 模块对 sentinel-1 SAR 数据进行预处理

一、前言

本章演示了用于 ESA SNAP 工具自动化的 Snappy Python 模块。

将展示用于预处理 Sentinel-1 SAR 数据的自动处理链的代码示例，包括 GRD（地面范围检测数据）的校准、子集和地形校正。

snappy的详细安装教程可以在这里找到：https://senbox.atlassian.net/wiki/display/SNAP/How+to+use+the+SNAP+API+from+Python

二、实现过程

（1）首先，导入所需的 Python 模块：

import snappy

from snappy import ProductIO
from snappy import HashMap

import os, gc   
from snappy import GPF

GPF.getDefaultInstance().getOperatorSpiRegistry().loadOperatorSpis()
HashMap = snappy.jpy.get_type('java.util.HashMap')

（2）现在遍历位于文件夹中的所有 Sentinel-1 数据子文件夹（当然，请确保数据已经解压缩）：

path = "D:\\SENTINEL\\"
 for folder in os.listdir(path):

   gc.enable()
   
   output = path + folder + "\\"  
   timestamp = folder.split("_")[4] 
   date = timestamp[:8]

（3）然后，读入 Sentinel-1 数据产品：

   sentinel_1 = ProductIO.readProduct(output + "\\manifest.safe")    
   print sentinel_1

如果偏振波段可用，请拆分您的代码以分别处理 VH 和 VV 强度数据。 第一步是将 DN 值分别转换为 Sigma Naught 的校准程序。 您也可以指定参数以分贝为单位输出图像。

pols = ['VH','VV'] 
   for p in pols:  
      polarization = p    
    
      ### CALIBRATION
  
      parameters = HashMap() 
      parameters.put('outputSigmaBand', True) 
      parameters.put('sourceBands', 'Intensity_' + polarization) 
      parameters.put('selectedPolarisations', polarization) 
      parameters.put('outputImageScaleInDb', False)  

      calib = output + date + "_calibrate_" + polarization 
      target_0 = GPF.createProduct("Calibration", parameters, sentinel_1) 
      ProductIO.writeProduct(target_0, calib, 'BEAM-DIMAP')

（4）接下来，指定子集 AOI 以减少数据量和处理时间。 AOI 由其外多边形指定，并使用 Well Known Text (WKT) 格式。

### SUBSET

      calibration = ProductIO.readProduct(calib + ".dim")    
      WKTReader = snappy.jpy.get_type('com.vividsolutions.jts.io.WKTReader')

      wkt = "POLYGON((12.76221 53.70951, 12.72085 54.07433, 13.58674 54.07981, 
                      13.59605 53.70875, 12.76221 53.70951))"

      geom = WKTReader().read(wkt)

      parameters = HashMap()
      parameters.put('geoRegion', geom)
      parameters.put('outputImageScaleInDb', False)

      subset = output + date + "_subset_" + polarization
      target_1 = GPF.createProduct("Subset", parameters, calibration)
      ProductIO.writeProduct(target_1, subset, 'BEAM-DIMAP')

（5）通过使用自动下载的 SRTM 3 arcsecond 产品 (90m)，应用距离多普勒地形校正来校正停留和透视缩短效应。 您还可以从本地路径指定具有更高空间分辨率的自己的 DEM 产品：

      ### TERRAIN CORRECTION
 
      parameters = HashMap()     
      parameters.put('demResamplingMethod', 'NEAREST_NEIGHBOUR') 
      parameters.put('imgResamplingMethod', 'NEAREST_NEIGHBOUR') 
      parameters.put('demName', 'SRTM 3Sec') 
      parameters.put('pixelSpacingInMeter', 10.0) 
      parameters.put('sourceBands', 'Sigma0_' + polarization)
 
      terrain = output + date + "_corrected_" + polarization 
      target_2 = GPF.createProduct("Terrain-Correction", parameters, subset) 
      ProductIO.writeProduct(target_2, terrain, 'GeoTIFF')