<jupyter_start><jupyter_text>ETHZ-01-03-01 Co-seismic deformation maps - Sentinel-1 As a data processor developer, I want to implement, and package an algorithm processing a pair of S1 SAR SLC datasets using the SNAP toolbox notebook archetype.  Quick link* [Objective](objective)* [Test Site](test-site)* [Context](context)* [Applicability](applicability)* [Data](data)* [Service Definition](service)* [Parameter Definition](parameter)* [Runtime Parameter Definition](runtime)* [Workflow](workflow)* [Strengths and Limitations](strengths-limitations) * [License](license)  Objective   Data SENTINEL data products are made available systematically and free of charge to all data users including the general public, scientific and commercial users. Radar data will be delivered within an hour of reception for Near Real-Time (NRT) emergency response, within three hours for NRT priority areas and within 24 hours for systematically archived data.All data products are distributed in the SENTINEL Standard Archive Format for Europe (SAFE) format.Data products are available in single polarisation (VV or HH) for Wave mode and dual polarisation (VV+VH or HH+HV) and single polarisation (HH or VV) for SM, IW and EW modes.  Service definition<jupyter_code>service = dict([('title', 'ETHZ-01-03-01 Co-seismic deformation maps'),
                ('abstract', 'This application takes a pair of Sentinel-1 products and generates and interferogram'),
                ('id', 'ewf-ethz-01-03-01')])<jupyter_output><empty_output><jupyter_text>Parameter Definition  Output file format:* BEAM-DIMAP* GeoTIFF-BigTIFF<jupyter_code>format = dict([('id', 'format'),
               ('value', 'GeoTIFF-BigTIFF'),
               ('title', 'Output file format'),
               ('abstract', 'Output file format: BEAM-DIMAP or GeoTIFF-BigTIFF')])<jupyter_output><empty_output><jupyter_text>Runtime parameter definition    The variable values in this section are only relevant for the basic test case. In an actual processing context, the values are replaced by those of the parameters for the process execution. **Input identifiers**This is the Sentinel-1 stack of master and slave products' identifiers:<jupyter_code>input_identifiers = ('S1A_IW_SLC__1SDV_20161030T163141_20161030T163208_013722_01603F_4094', 'S1A_IW_SLC__1SDV_20161018T163141_20161018T163208_013547_015AEB_5994')<jupyter_output><empty_output><jupyter_text>**Input references**This is the Sentinel-1 stack of catalogue references:<jupyter_code>input_references = ('https://catalog.terradue.com/sentinel1/search?format=atom&uid=S1A_IW_SLC__1SDV_20161030T163141_20161030T163208_013722_01603F_4094', 'https://catalog.terradue.com/sentinel1/search?format=atom&uid=S1A_IW_SLC__1SDV_20161018T163141_20161018T163208_013547_015AEB_5994')<jupyter_output><empty_output><jupyter_text>**Data path**This path defines where the data is staged-in:<jupyter_code>data_path = "/workspace/data"<jupyter_output><empty_output><jupyter_text>Workflow  Import the packages required for processing the data<jupyter_code>from snappy import jpy
from snappy import ProductIO
from snappy import GPF
from snappy import HashMap

import dateutil.parser as parser
import gc
from datetime import datetime

import cioppy

import gdal
import osr

from shapely.wkt import loads
from shapely.geometry import box

import lxml.etree as etree
import numpy as np
import os

import shutil
import lxml.etree as etree
import subprocess
import tempfile
import time
import psutil
from snappy import jpy
from snappy import ProductIO
from snappy import GPF
from snappy import HashMap
import os

class GraphProcessor():
    
    def __init__(self):
        self.root = etree.Element('graph')
    
        version = etree.SubElement(self.root, 'version')
        version.text = '1.0'
        self.pid = None
        self.p = None
   
    def view_graph(self):
        
        print etree.tostring(self.root , pretty_print=True)
        
    def add_node(self, node_id, operator, parameters, source):
    
        xpath_expr = '/graph/node[@id="%s"]' % node_id

        if len(self.root.xpath(xpath_expr)) != 0:

            node_elem = self.root.xpath(xpath_expr)[0]
            operator_elem = self.root.xpath(xpath_expr + '/operator')[0]
            sources_elem = self.root.xpath(xpath_expr + '/sources')[0]
            parameters_elem = self.root.xpath(xpath_expr + '/parameters')

            for key, value in parameters.iteritems():
                p_elem = self.root.xpath(xpath_expr + '/parameters/%s' % key)[0]
                p_elem.text = value
        else:

            node_elem = etree.SubElement(self.root, 'node')
            operator_elem = etree.SubElement(node_elem, 'operator')
            sources_elem = etree.SubElement(node_elem, 'sources')

            if isinstance(source, list):

                for index, s in enumerate(source):
                    if index == 0:  
                        source_product_elem = etree.SubElement(sources_elem, 'sourceProduct')

                    else: 
                        source_product_elem = etree.SubElement(sources_elem, 'sourceProduct.%s' % str(index))

                    source_product_elem.attrib['refid'] = s

            elif source != '':
                source_product_elem = etree.SubElement(sources_elem, 'sourceProduct')
                source_product_elem.attrib['refid'] = source

            parameters_elem = etree.SubElement(node_elem, 'parameters')
            parameters_elem.attrib['class'] = 'com.bc.ceres.binding.dom.XppDomElement'

            for key, value in parameters.iteritems():

                parameter_elem = etree.SubElement(parameters_elem, key)
                parameter_elem.text = value

        node_elem.attrib['id'] = node_id

        operator_elem.text = operator 

    def save_graph(self, filename):
        
        with open(filename, 'wb') as file:
            file.write('<?xml version="1.0" encoding="UTF-8"?>\n')
            file.write(etree.tostring(self.root, pretty_print=True))
     
    def plot_graph(self):
        
        for node_id in self.root.xpath('/graph/node/@id'):
            

            xpath_expr = '/graph/node[@id="%s"]' % node_id
            
            if len(self.root.xpath(xpath_expr + '/sources/sourceProduct')) != 0:
                print(self.root.xpath(xpath_expr + '/sources/sourceProduct'))[0].attrib['refid']
                print node_id
            else:
                print node_id
        return True
    
    def run(self):
        
        fd, path = tempfile.mkstemp()
        
        try:
        
            self.save_graph(filename=path)
            options = ['/opt/snap6/bin/gpt',
               '-x',
               '-c',
               '2048M',
               path]

            p = subprocess.Popen(options,
                stdout=subprocess.PIPE, stdin=subprocess.PIPE, stderr=subprocess.PIPE)

            print p.pid
            res, err = p.communicate()
            print res, err
            if p.returncode != 0:
                raise Exception('An error occurred during the execution of gpt (see log)')
            
        except Exception as e:
            with open('stdout.txt', 'wb') as file:
                file.write(res)
                file.close()
            with open('stderr.txt', 'wb') as file:
                file.write(err)
                file.close()
            
            raise
        finally:
            os.remove(path)
        
def get_snap_parameters(operator):
    
    op_spi = GPF.getDefaultInstance().getOperatorSpiRegistry().getOperatorSpi(operator)

    op_params = op_spi.getOperatorDescriptor().getParameterDescriptors()

    return op_params
mygraph = GraphProcessor()<jupyter_output><empty_output><jupyter_text>Read the products<jupyter_code>operator = 'Read'

read_nodes = []

for index, identifier in enumerate(input_identifiers):
    
    parameters = dict()

    for param in get_snap_parameters(operator):
    
        if param.getName() == 'file':
            parameters[param.getName()] = os.path.join(data_path, identifier, identifier + '.SAFE', 'manifest.safe')    
        else:
            parameters[param.getName()] = param.getDefaultValue()
    node_id = 'Read(%s)' % index
    
    read_nodes.append(node_id)
    
    mygraph.add_node(node_id, 'Read', parameters, '')<jupyter_output><empty_output><jupyter_text>Apply orbit file<jupyter_code>operator = 'Apply-Orbit-File'

parameters = dict()

for param in get_snap_parameters(operator):
    
    print(param.getName(), param.getDefaultValue())
    
    parameters[param.getName()] = param.getDefaultValue()
orbit_nodes = []

for index, source_node in enumerate(read_nodes):
    
    node_id = 'Apply-Orbit-File(%s)' % index
    
    orbit_nodes.append(node_id)
    
    mygraph.add_node(node_id, 'Apply-Orbit-File', parameters, source_node)<jupyter_output><empty_output><jupyter_text>TOPSAR split<jupyter_code>operator = 'TOPSAR-Split'

parameters = dict()

for param in get_snap_parameters(operator):
    
    print(param.getName(), param.getDefaultValue())
    
    parameters[param.getName()] = param.getDefaultValue()
slave_split_nodes = []
master_split_nodes = []

for index_node, source_node in enumerate(orbit_nodes):
    
    for index, subswath in enumerate(['IW1', 'IW2', 'IW3']):  
    
        parameters['subswath'] =  subswath
        parameters['selectedPolarisations'] = 'VV'

        node_id = 'TOPSAR-Split(%s)' % str(index_node * 3 + index)
    
        if index_node == 0:
            slave_split_nodes.append(node_id)
        else:
            master_split_nodes.append(node_id)
    
        mygraph.add_node(node_id, operator, parameters, source_node)<jupyter_output><empty_output><jupyter_text>Back-Geocoding<jupyter_code>operator = 'Back-Geocoding'

parameters = dict()

for param in get_snap_parameters(operator):
    
    print(param.getName(), param.getDefaultValue())
    
    parameters[param.getName()] = param.getDefaultValue()
backgeo_nodes = []
    
for index, subswath in enumerate(['IW1', 'IW2', 'IW3']):  
            
    node_id = 'Back-Geocoding(%s)' % index
   
    source_nodes = [slave_split_nodes[index], master_split_nodes[index]]
    
    mygraph.add_node(node_id, operator, parameters, source_nodes)

    backgeo_nodes.append(node_id)<jupyter_output><empty_output><jupyter_text>Interferogram<jupyter_code>operator = 'Interferogram'

parameters = dict()

for param in get_snap_parameters(operator):
    
    print(param.getName(), param.getDefaultValue())
    
    parameters[param.getName()] = param.getDefaultValue()
interferogram_nodes = []

for index, subswath in enumerate(['IW1', 'IW2', 'IW3']):  
            
    node_id = 'Interferogram(%s)' % index
   
    source_node = backgeo_nodes[index]
    
    mygraph.add_node(node_id, operator, parameters, source_node)

    interferogram_nodes.append(node_id)<jupyter_output><empty_output><jupyter_text>TOPSAR Deburst<jupyter_code>operator = 'TOPSAR-Deburst'

parameters = dict()

for param in get_snap_parameters(operator):
    
    print(param.getName(), param.getDefaultValue())
    
    parameters[param.getName()] = param.getDefaultValue()
deburst_nodes = []

for index, subswath in enumerate(['IW1', 'IW2', 'IW3']):  
 
    parameters['selectedPolarisations'], 'VV'

    node_id = 'TOPSAR-Deburst(%s)' % index
   
    source_node = interferogram_nodes[index]
    
    mygraph.add_node(node_id, operator, parameters, source_node)

    deburst_nodes.append(node_id)<jupyter_output><empty_output><jupyter_text>TOPSAR Merge<jupyter_code>operator = 'TOPSAR-Merge'

parameters = dict()

for param in get_snap_parameters(operator):
    
    print(param.getName(), param.getDefaultValue())
    
    parameters[param.getName()] = param.getDefaultValue()
parameters['selectedPolarisations'], 'VV'

source_nodes = deburst_nodes

node_id = operator

mygraph.add_node(node_id, operator, parameters, source_nodes)<jupyter_output><empty_output><jupyter_text>TopoPhaseRemoval<jupyter_code>operator = 'TopoPhaseRemoval'

parameters = dict()

for param in get_snap_parameters(operator):
    
    print(param.getName(), param.getDefaultValue())
    
    parameters[param.getName()] = param.getDefaultValue()
source_node = 'TOPSAR-Merge'

node_id = operator

mygraph.add_node(node_id, operator, parameters, source_node)<jupyter_output><empty_output><jupyter_text>GoldsteinPhaseFiltering<jupyter_code>operator = 'GoldsteinPhaseFiltering'

parameters = dict()

for param in get_snap_parameters(operator):
    
    print(param.getName(), param.getDefaultValue())
    
    parameters[param.getName()] = param.getDefaultValue()
source_node = 'TopoPhaseRemoval'

node_id = operator

mygraph.add_node(node_id, operator, parameters, source_node)<jupyter_output><empty_output><jupyter_text>Terrain correction<jupyter_code>operator = 'Terrain-Correction'

parameters = dict()

for param in get_snap_parameters(operator):
    
    print(param.getName(), param.getDefaultValue())
    
    parameters[param.getName()] = param.getDefaultValue()
source_node = 'GoldsteinPhaseFiltering'

parameters['outputComplex'] = 'true'

node_id = operator

mygraph.add_node(node_id, operator, parameters, source_node)
ciop = cioppy.Cioppy()

search = ciop.search(end_point=input_references[0],
                     params=[],
                     output_fields='enclosure,identifier,startdate,enddate,wkt,orbitNumber,orbitDirection,swathIdentifier,wrsLongitudeGrid',
                     model='EOP') 

search2 = ciop.search(end_point=input_references[1],
                     params=[],
                     output_fields='enclosure,identifier,startdate,enddate,wkt,orbitNumber,orbitDirection,swathIdentifier,wrsLongitudeGrid',
                     model='EOP')

if (search[0]['startdate'] < search2[0]['startdate']):
    start_date = search[0]['startdate']
else:
    start_date = search2[0]['startdate']

if (search[0]['enddate'] > search2[0]['enddate']):
    end_date = search[0]['enddate']
else:
    end_date = search2[0]['enddate']
output_name = 'S1_IFG_%s_%s_%s' % ("%03d"%int(search[0]['wrsLongitudeGrid']), 
                                    parser.parse(start_date).strftime('%Y%m%d_%H%M%S'),
                                    parser.parse(end_date).strftime('%Y%m%d_%H%M%S'))
operator = 'Write'
 
parameters = dict()

for param in get_snap_parameters(operator):
    
    if param.getName() == 'file':
        
        param_value = output_name
             
    elif param.getName() == 'formatName':
                
        param_value = 'BEAM-DIMAP'
        
    else:
    
        param_value = param.getDefaultValue()
    
    
    print param.getName(), param_value
    
    parameters[param.getName()] = param_value


mygraph.add_node(operator, 
             operator, 
             parameters,
             'Terrain-Correction')
mygraph.view_graph()
mygraph.run()
reader = ProductIO.getProductReader("BEAM-DIMAP")

ifg = reader.readProductNodes(output_name + '.dim', None)
band_names = ifg.getBandNames()
target_bands = list(band_names) 
target_bands.append(target_bands[0].replace('i_', 'phase_'))
BandDescriptor = jpy.get_type('org.esa.snap.core.gpf.common.BandMathsOp$BandDescriptor')

targetBands = jpy.array('org.esa.snap.core.gpf.common.BandMathsOp$BandDescriptor', len(target_bands))

for index, band in enumerate(target_bands):
     
    targetBand = BandDescriptor()
    
    if 'phase_' in band:
       
        targetBand.expression = 'atan2(%s, %s)' % (target_bands[1], target_bands[0])
   
    else:
    
        targetBand.expression = band

    targetBand.name = band
    targetBand.type = 'float32'
    
    targetBands[index]= targetBand
        

parameters = HashMap()
parameters.put('targetBands', targetBands)

result = GPF.createProduct('BandMaths', parameters, ifg)
ProductIO.writeProduct(result, 
                       output_name + '.tif',
                       'GeoTIFF-BigTiff')
def eop_metadata(metadata):

    opt = 'http://www.opengis.net/opt/2.1'
    om  = 'http://www.opengis.net/om/2.0'
    gml = 'http://www.opengis.net/gml/3.2'
    eop = 'http://www.opengis.net/eop/2.1'
    sar = 'http://www.opengis.net/sar/2.1'
    
    root = etree.Element('{%s}EarthObservation' % sar)

    phenomenon_time = etree.SubElement(root, '{%s}phenomenonTime' % om)

    time_period = etree.SubElement(phenomenon_time, '{%s}TimePeriod' % gml)

    begin_position = etree.SubElement(time_period, '{%s}beginPosition'  % gml)

    end_position = etree.SubElement(time_period, '{%s}endPosition'  % gml)

    procedure = etree.SubElement(root, '{%s}procedure' % om)

    earth_observation_equipment = etree.SubElement(procedure, '{%s}EarthObservationEquipment' % eop)

    acquisition_parameters = etree.SubElement(earth_observation_equipment, '{%s}acquisitionParameters' % eop)

    acquisition = etree.SubElement(acquisition_parameters, '{%s}Acquisition' % sar)

    orbit_number = etree.SubElement(acquisition, '{%s}orbitNumber' % eop)

    wrs_longitude_grid = etree.SubElement(acquisition, '{%s}wrsLongitudeGrid' % eop)

    polarisation_channels = etree.SubElement(acquisition, '{%s}polarisationChannels' % eop)
    
    feature_of_interest = etree.SubElement(root, '{%s}featureOfInterest' % om)
    footprint = etree.SubElement(feature_of_interest, '{%s}Footprint' % eop)
    multi_extentOf = etree.SubElement(footprint, '{%s}multiExtentOf' % eop)
    multi_surface = etree.SubElement(multi_extentOf, '{%s}MultiSurface' % gml)
    surface_members = etree.SubElement(multi_surface, '{%s}surfaceMembers' % gml)
    polygon = etree.SubElement(surface_members, '{%s}Polygon' % gml)    
    exterior = etree.SubElement(polygon, '{%s}exterior' % gml)  
    linear_ring = etree.SubElement(exterior, '{%s}LinearRing' % gml) 
    poslist = etree.SubElement(linear_ring, '{%s}posList' % gml) 


    result = etree.SubElement(root, '{%s}result' % om)
    earth_observation_result = etree.SubElement(result, '{%s}EarthObservationResult' % opt)
    cloud_cover_percentage = etree.SubElement(earth_observation_result, '{%s}cloudCoverPercentage' % opt)
    
    metadata_property = etree.SubElement(root, '{%s}metaDataProperty' % eop)
    earth_observation_metadata = etree.SubElement(metadata_property, '{%s}EarthObservationMetaData' % eop)
    identifier = etree.SubElement(earth_observation_metadata, '{%s}identifier' % eop)
    
    begin_position.text = metadata['startdate']
    end_position.text = metadata['enddate']
    #orbit_number.text = metadata['orbitNumber']
    wrs_longitude_grid.text = metadata['wrsLongitudeGrid']
    
    coords = np.asarray([t[::-1] for t in list(loads(metadata['wkt']).exterior.coords)]).tolist()
 
    pos_list = ''
    for elem in coords:
        pos_list += ' '.join(str(e) for e in elem) + ' '   

    poslist.attrib['count'] = str(len(coords))
    poslist.text = pos_list
    
    
    identifier.text = metadata['identifier'] 

    return etree.tostring(root, pretty_print=True)
src = gdal.Open(output_name + '.tif')
ulx, xres, xskew, uly, yskew, yres  = src.GetGeoTransform()

max_x = ulx + (src.RasterXSize * xres)
min_y = uly + (src.RasterYSize * yres)
min_x = ulx 
max_y = uly

source = osr.SpatialReference()
source.ImportFromWkt(src.GetProjection())

target = osr.SpatialReference()
target.ImportFromEPSG(4326)

transform = osr.CoordinateTransformation(source, target)

result_wkt = box(transform.TransformPoint(min_x, min_y)[0],
        transform.TransformPoint(min_x, min_y)[1],
        transform.TransformPoint(max_x, max_y)[0],
        transform.TransformPoint(max_x, max_y)[1]).wkt
search[0]['identifier'] = output_name
search[0]['wkt'] = result_wkt

search[0]['startdate'] = start_date
search[0]['enddate'] = end_date

eop_xml = output_name + '.xml'
with open(eop_xml, 'wb') as file:
    file.write('<?xml version="1.0" encoding="UTF-8"?>\n')
    file.write(eop_metadata(search[0]))
for properties_file in ['result', 'stage-in']:

    if properties_file == 'result':
        title = 'Reproducibility notebook used for generating %s' % output_name
    else: 
        title = 'Reproducibility stage-in notebook for Sentinel-1 data for generating %s' % output_name
        
    with open(properties_file + '.properties', 'wb') as file:
        file.write('title=%s\n' % title)
        file.write('date=%s/%s\n' % (start_date, end_date))
        file.write('geometry=%s' % (search[0]['wkt']))<jupyter_output><empty_output><jupyter_text>Clean-up<jupyter_code>shutil.rmtree(output_name + '.data')
os.remove(output_name + '.dim')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12099609375, 'fraction_numerical': 0.01416015625, 'mean_word_length': 3.5291906236178683, 'pattern_counts': {'":': 0, '<': 72, '<?xml version=': 2, '>': 72, 'https://': 2, 'lorem ipsum': 0, 'www.': 5, 'xml': 7}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '168409', 'license': 'Apache-2.0', 'n_tokens_mistral': 6826, 'n_tokens_neox': 6552, 'n_words': 1224, 'path': 'src/main/app-resources/notebook/libexec/input.ipynb', 'repo_name': 'ec-better/ewf-ethz-01-03-01', 'size': '20480'}