Creating Mosaic of Sentinel-1 using Snappy

In this tutorial, we are going to perform the mosaicing of two adjacent Sentinel-1 scenes using Snappy, the Python interface for SNAP.

If you haven’t, follow along this tutorial to see how you can setup the development environment for SNAP in your machine. To overview, the basic steps that we cover in this tutorial are:

1. Apply Orbit File
2. Border Noise Removal from scenes
3.  Calibration of each scene
4. Multilooking
5. Range Doppler Terrain Correction using DEM
6. Mosaicing of two or more scenes

Use the script below to perform mosaicing:

import snappy
from snappy import ProductIO, GPF, HashMap
import os
from functools import partial

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

input = 'data'
output = 'output'

def write_product(data, file_path, format=None):
    # allowed format are 
    # GeoTIFF-BigTIFF,HDF5,Snaphu,BEAM-DIMAP,GeoTIFF+XML,PolSARPro,NetCDF-CF,NetCDF-BEAM,ENVI,JP2,Generic Binary BSQ,Gamma,CSV,NetCDF4-CF,GeoTIFF,NetCDF4-BEAM
    ProductIO.writeProduct(data, file_path, format if format else 'BEAM-DIMAP')

def apply_orbit_file(data, datestamp):
    # uncomment and indent if you want to write the output
    #if not os.path.isfile(os.path.join(output, '{}_Orb.dim').format(datestamp)) and not os.path.isfile(os.path.join(output, '{}_Orb.data').format(datestamp)):
    params = HashMap()
    # find a way to put this in except
    params.put('continueOnFail', False)
    orbit = GPF.createProduct('Apply-Orbit-File', params, data)
    #write_product(orbit, os.path.join(output, '{}_Orb'.format(datestamp)))
    return orbit

def border_noise_removal(data, datestamp):
    params = HashMap()
    # put 1000 to make sure every bad pixels are included
    params.put('borderLimit', 1000)
    params.put('trimThreshold', 0.5)
    noise_removed = GPF.createProduct('Remove-GRD-Border-Noise', params, data)
    #write_product(noise_removed, os.path.join(output, '{}_Orb_Bdr'.format(datestamp)))
    return noise_removed

def do_calibration(noise_removed, datestamp):
    '''
    orbit: orbital corrected file
    '''

    params = HashMap()
    params.put('outputSigmaBand', True)
    params.put('outputGammaBand', False)
    params.put('outputBetaBand', False)
    calibration = GPF.createProduct('Calibration', params, noise_removed)
    #write_product(calibration, os.path.join(output, '{}_Orb_Bdr_Cal'.format(datestamp)))
    return calibration

def perform_multilook(calibration, datestamp, range_look_number=3, azimuth_look_number=3):
    params = HashMap()
    params.put('nRgLooks', range_look_number)
    params.put('nAzLooks', azimuth_look_number)
    params.put('outputIntensity', True)
    params.put('grSquarePixel', True)
    multilook = GPF.createProduct('Multilook', params, calibration)
    #write_product(multilook, os.path.join(output, '{}_Orb_Bdr_Cal_ML'.format(datestamp)))
    return multilook

def range_doppler_terrain_correction(multilook, datestamp):
    params = HashMap()
    params.put('demName', 'SRTM 1Sec HGT')
    params.put('pixelSpacingInMeter', 30.0)
    params.put('pixelSpacingInDegree', 0.00026949458)
    terrain = GPF.createProduct('Terrain-Correction', params, multilook)
    write_product(terrain, os.path.join(output, '{}_Orb_Bdr_Cal_ML_TC'.format(datestamp)))
    return terrain

def sar_mosaic(terrains):
    print('started mosaic')
    params = HashMap()
    params.put('pixelSize', 30.0)
    mosaic = GPF.createProduct('SAR-Mosaic', params, terrains)
    write_product(mosaic, os.path.join(output, 'Mosaic'))
    return mosaic

def process(file):
    data = ProductIO.readProduct(os.path.join(input, file))
    datestamp = file.split('_')[4][:8]
    orbit = apply_orbit_file(data, datestamp)
    noise_removed = border_noise_removal(data, datestamp)
    calibration = do_calibration(noise_removed, datestamp)
    multilook = perform_multilook(calibration, datestamp)
    terrain = range_doppler_terrain_correction(multilook, datestamp)
    return terrain

def set_path():
    os.path.dirname(os.path.dirname(__file__))
    path = os.path.join(os.getcwd())
    os.chdir(path)
    return path

def main():
    set_path()
    files = [f for f in os.listdir(input) if f.endswith('.zip')]
    params = snappy.jpy.get_type('java.util.HashMap')
    terrain_corrected_files = []
    for file in files:
        print('started file: {}'.format(file))
        terrain = process(file)
        terrain_corrected_files.append(terrain)
    mosaic = sar_mosaic(terrain_corrected_files)
    print('finished mosaicing')

if __name__ == '__main__':
    main()