Step 4 - Implement the data transformation steps

The data transformation application will take one Sentinel-1 products (staged-in during Step 3 - Stage-in the EO data of this tutorial) and generate the backscatter using the Sentinel Application Platform (SNAP) Python bindings.

The data transformation steps above will use a Jupyter Notebook to:

• Read the Sentinel-1 product using the SNAP reader classes
• Application of orbit file using the SNAP Apply-Orbit-File operator
• Border noise removal using the SNAP ThermalNoiseRemoval operator
• Calibration using the SNAP Calibration operator
• Speckle filtering using the SNAP Speckle-Filter operator
• Terrain correction using the SNAP Terrain-Correction operator
• Conversion to dB using the SNAP linearToFromdB operator

At runtime this Jupyter Notebook will be instantiated and parametrized to process the Sentinel-1 archive.

Service definition

The data transformation application will be exposed as a Web Processing Service.

The definition of the Web Processing Service information such as the title and the abstract is done with a Python dictionary:

service = dict([('title', 'WFP-01-01-01 Sentinel-1 backscatter timeseries'),
                ('abstract', 'WFP-01-01-01 Data transformation application - Sentinel-1 backscatter timeseries'),
                ('id', 'wfp-01-01-01')])

Note: As you will see in Step 8 - Deploy the application on the Production Centre, the Web Processing Service process identifier is set by the Production Centre. The ‘id’ set above is used to provide and identifier to the workflow

Parameter Definition

Data transformation application may have to expose parameters that can be changed via the Web Processing Service interface at submission time.

These parameters are defined using a Python dictionary that defines the parameter identifier, its title and abstract and finally its default value:

filterSizeX = dict([('id', 'filterSizeX'),
               ('value', '5'),
               ('title', 'Speckle-Filter filterSizeX'),
               ('abstract', 'Set the Speckle-Filter filterSizeX (defaults to 5)')])

To use the parameter value, symply do:

int(filterSizeX['value'])

5

Runtime parameter definition

Runtime parameters are mandatory and define those parameters whose values will be changed at runtime.

These are:

• input_identifier - this is the Sentinel-1 product identifier. At runtime its value is replaced with the Sentinel-1 product identifier being processed
• input_reference - this is the Sentinel-1 product catalogue entry URL. At runtime its value is also replaced with the Sentinel-1 product catalogue entry URL being processed
• data_path - this is the local path where the Sentinel-1 was staged-in in Step 3 - Stage-in the EO data. At runtime its value is replaced by a folder with an unique value

Publishing outputs or results

To publish results simply save the output or the result locally.

Implementing the data transformation processing steps

Note: At this stage, you can either copy each of the code cells to the input.ipynb or download the complete input.ipynb with the URL provided at the bottom of this page

• Define the Service definition with:

service = dict([('title', 'WFP-01-01-01 Sentinel-1 backscatter timeseries'),
                ('abstract', 'WFP-01-01-01 Data transformation application - Sentinel-1 backscatter timeseries'),
                ('id', 'wfp-01-01-01')])

• Define the Web Processing Service parameters:

filterSizeX = dict([('id', 'filterSizeX'),
               ('value', '5'),
               ('title', 'Speckle-Filter filterSizeX'),
               ('abstract', 'Set the Speckle-Filter filterSizeX (defaults to 5)')])

filterSizeY = dict([('id', 'filterSizeY'),
               ('value', '5'),
               ('title', 'Speckle-Filter filterSizeY'),
               ('abstract', 'Set the Speckle-Filter filterSizeY (defaults to 5)')])

polarisation = dict([('id', 'polarisation'),
               ('value', 'VV'),
               ('title', 'Sentinel-1 polarisation (VV or HH)'),
               ('abstract', 'Sentinel-1 polarisation (VV or HH)')])

wkt = dict([('id', 'wkt'),
            ('value', 'POLYGON((-5.5 17.26, -1.08 17.26, -1.08 13.5, -5.5 13.5, -5.5 17.26))'),
            ('title', 'Area of interest in WKT'),
            ('abstract', 'Area of interest using a polygon in Well-Known-Text format')])

• Define the Runtime parameters

input_identifier = 'S1A_IW_GRDH_1SDV_20171210T182024_20171210T182049_019644_021603_0A33'

input_reference = 'https://catalog.terradue.com/sentinel1/search?format=atom&uid=%s' % input_identifier

data_path = '/workspace/data'

• Define the Python libraries

%matplotlib inline

import warnings
warnings.filterwarnings("ignore")
import os
import sys
import glob
sys.path.append('/opt/anaconda/bin/')

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as colors

from snappy import jpy
from snappy import ProductIO
from snappy import GPF
from snappy import HashMap

import gc

• Read the Sentinel-1 product

s1meta = "manifest.safe"

s1prd = os.path.join(data_path, input_identifier, input_identifier + '.SAFE', s1meta)

reader = ProductIO.getProductReader("SENTINEL-1")
product = reader.readProductNodes(s1prd, None)

• Apply the Thermal Noise Removal

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

parameters = HashMap()

parameters.put('selectedPolarisations', polarisation['value'])
parameters.put('removeThermalNoise', 'true')
parameters.put('reIntroduceThermalNoise', 'false')

thermal_noise_removal = GPF.createProduct('ThermalNoiseRemoval', parameters, product)

Note: You might be wondering how to discover what parameters should be used with a given SNAP operator. The cell below shows how to do it programmatically

operator = 'ThermalNoiseRemoval'

op_spi = GPF.getDefaultInstance().getOperatorSpiRegistry().getOperatorSpi(operator)

op_params = op_spi.getOperatorDescriptor().getParameterDescriptors()

for param in op_params:
    print(param.getName(), param.getDefaultValue())

('selectedPolarisations', None)
('removeThermalNoise', 'true')
('reIntroduceThermalNoise', 'false')

• Apply the orbit file correction

parameters = HashMap()

parameters.put('orbitType', 'Sentinel Precise (Auto Download)')
parameters.put('polyDegree', '3')
parameters.put('continueOnFail', 'false')

apply_orbit_file = GPF.createProduct('Apply-Orbit-File', parameters, thermal_noise_removal)

• Perform the Calibration

parameters = HashMap()

parameters.put('auxFile', 'Product Auxiliary File')
parameters.put('outputImageInComplex', 'false')
parameters.put('outputImageScaleInDb', 'false')
parameters.put('createGammaBand', 'false')
parameters.put('createBetaBand', 'true')
parameters.put('selectedPolarisations', polarisation['value'])
parameters.put('outputSigmaBand', 'false')
parameters.put('outputGammaBand', 'false')
parameters.put('outputBetaBand', 'true')

calibration = GPF.createProduct('Calibration', parameters, apply_orbit_file)

• Apply the Speckle Filter

parameters = HashMap()

parameters.put('sourceBands', 'Beta0_%s' % (polarisation['value']))
parameters.put('filter', 'Lee')
parameters.put('filterSizeX', filterSizeX['value'])
parameters.put('filterSizeY', filterSizeY['value'])
parameters.put('dampingFactor', '2')
parameters.put('estimateENL', 'true')
parameters.put('enl', '1.0')
parameters.put('numLooksStr', '1')
parameters.put('targetWindowSizeStr', '3x3')
parameters.put('sigmaStr', '0.9')
parameters.put('anSize', '50')

speckle_filter = GPF.createProduct('Speckle-Filter', parameters, calibration)

• Apply the Terrain Correction

parameters = HashMap()

parameters.put('sourceBands', 'Beta0_%s' % (polarisation['value']))
parameters.put('demName', 'SRTM 3Sec')
parameters.put('externalDEMFile', '')
parameters.put('externalDEMNoDataValue', '0.0')
parameters.put('externalDEMApplyEGM', 'true')
parameters.put('demResamplingMethod', 'BILINEAR_INTERPOLATION')
parameters.put('imgResamplingMethod', 'BILINEAR_INTERPOLATION')
parameters.put('pixelSpacingInMeter', '10.0')
#parameters.put('pixelSpacingInDegree', '8.983152841195215E-5')
parameters.put('mapProjection', 'AUTO:42001')
parameters.put('nodataValueAtSea', 'true')
parameters.put('saveDEM', 'false')
parameters.put('saveLatLon', 'false')
parameters.put('saveIncidenceAngleFromEllipsoid', 'false')
parameters.put('saveProjectedLocalIncidenceAngle', 'false')
parameters.put('saveSelectedSourceBand', 'true')
parameters.put('outputComplex', 'false')
parameters.put('applyRadiometricNormalization', 'false')
parameters.put('saveSigmaNought', 'false')
parameters.put('saveGammaNought', 'false')
parameters.put('saveBetaNought', 'false')
parameters.put('incidenceAngleForSigma0', 'Use projected local incidence angle from DEM')
parameters.put('incidenceAngleForGamma0', 'Use projected local incidence angle from DEM')
parameters.put('auxFile', 'Latest Auxiliary File')

terrain_correction = GPF.createProduct('Terrain-Correction', parameters, speckle_filter)

• Apply the Linar to dB

parameters = HashMap()

lineartodb = GPF.createProduct('linearToFromdB', parameters, terrain_correction)

• Free some memory with the garbage collector

thermal_noise_removal = None
apply_orbit_file = None
calibration = None
speckle_filter = None
terrain_correction = None

gc.collect()

6

• Save the result

output_name = '%s_Beta0_%s.tif' % (input_identifier, polarisation['value'])

ProductIO.writeProduct(lineartodb,
                       output_name,
                       'GeoTIFF-BigTIFF')

Note: The data transformation is done!

Now that we have implemented the data transformation, let’s plot a subset of the backscatter:

list(lineartodb.getBandNames())

parameters = HashMap()

parameters.put('sourceBands', 'Beta0_%s_db' % (polarisation['value']))

parameters.put('geoRegion', 'POLYGON((-4.5 16.26, -2.08 16.26, -2.08 14.5, -4.5 14.5, -4.5 16.26))')
parameters.put('subSamplingX', '1')
parameters.put('subSamplingY', '1')
parameters.put('fullSwath', 'false')
parameters.put('tiePointGridNames', '')
parameters.put('copyMetadata', 'true')

subset = GPF.createProduct('Subset',
                           parameters,
                           lineartodb)

list(subset.getBandNames())

def plotBand(product, band, vmin, vmax):

    band = product.getBand(band)

    w = band.getRasterWidth()
    h = band.getRasterHeight()

    band_data = np.zeros(w * h, np.float32)
    band.readPixels(0, 0, w, h, band_data)

    band_data.shape = h, w

    imgplot = plt.imshow(band_data,
                         cmap=plt.cm.binary_r,
                         vmin=vmin,
                         vmax=vmax)


    return imgplot

plotBand(subset,
         list(lineartodb.getBandNames())[0],
         0,
         1000)

Documenting the Jupyter Notebook streaming executable

One of the nice features of Jupyter Notebooks is that they can incorporate text that documents what is done.

Download this link for a proposal for documenting the WFP-01-01-01 Sentinel-1 backscatter timeseries data transformation streaming notebook.

You can use it as the notebook streaming executable.

To do so, use the Jupyter left panel to browse to /workspace/wfp-01-01-01/src/main/app-resources/notebook/libexec and upload the ‘input.ipynb’ file after extracting it.

The next step will deploy the data transformation as a local Web Processing Servic and test it against a Sentinel-1 product using the Sandbox resources