Aggregation and Interpolation of Sentinel products - should I use snappy or GDAL tools?

You can open the .jp2 files with GDAL. At least L1A or the ones produced with sen2cor 2.0.6. . When trying gdal on L2A .jp2 from sen2cor 2.2.1 i get an error

To aggregate data you can use for example the Resampling operator or you can use the Binning operator.
You can invoke them like
result = GPF.createProduct('Resampling', parameters, product)
or
result = GPF.createProduct('Binning', parameters, products)
In general I think you can do the work with SNAP. You can try your workflow in the Desktop application first and if it works for you can also do it from snappy. I’ve updated the post http://forum.step.esa.int/t/about-the-python-category to contain links to some documentation. Also we will be glad to support you via this forum if you have specific problems .

To give you more hints it would be good to know your use case in more detail.

Additionally, it might be sufficient for you if you just configure some operators and simply run a graph with gpt.

@marpet Thank you for answer, but I have problem with the following part:

You can try your workflow in the Desktop application first and if it works for you can also do it from snappy

When I use SNAP desktop application, I don’t see any code (i.e. no Python/snappy code) in such a way that for example GRASS shows. How can I know how to execute exactly the same tool using snappy/Python?

@marpet I discovered that ‘full name’ of the SNAP plugin/operator which I need to use in Python code is within the help (as I mentioning in http://forum.step.esa.int/t/obtaining-sigma-backscatter-from-sentinel-1-subset.

I still have issues, but at least I know how to start.

You don’t need to use the full name.
The name should be sufficient if you use GPF as in this example

GPF.createProduct('BandMaths', parameters, product)
BandMaths is here the name of the operator.

In the help you can find a list of some operators, the basic ones.
But if you go on the command line you can type gpt -h and get the full list of all available operators.
I think in most cases it is possible to guess which processor from the Desktop application maps to which operator. But by looking up the name in the about box you are sure.
Except for the subset operation, this you can’t look up. But the name is simply ‘Subset

Sorry, I should have answered before.

@marpet Great, but this way requires to prepare proper “parameters” before calling an operator.
To be honest, I don’t know what are the available/needed (and how to find them) list or arrays of parameters for ‘Resampling’ or ‘Binning’.

When I typed full name of the operator in search engine, I proceeded to the documentation:
http://step.esa.int/docs/v3.0/apidoc/engine/org/esa/snap/core/gpf/common/resample/ResamplingOp.html
but I can read here only that ‘setParameter’ method is inherited from Operator class.
Still no idea how to set parameters for Resampling.

Hi!

I do this:

[me@sb-10-16-10-54 ~]$ /opt/snap-3.0/bin/gpt -h Terrain-Correction

which dumps the help information for the Terrain-Correction operator (just an example), then I use that info to create the parameters hashmap and then invoke the operator:

HashMap = jpy.get_type('java.util.HashMap') GPF.getDefaultInstance().getOperatorSpiRegistry().loadOperatorSpis() parameters = HashMap() parameters.put('demName', 'SRTM 3Sec') parameters.put('externalDEMNoDataValue', 0.0) parameters.put('demResamplingMethod', "BILINEAR_INTERPOLATION") parameters.put('imgResamplingMethod', "BILINEAR_INTERPOLATION") parameters.put('pixelSpacingInMeter', 10.0) parameters.put('pixelSpacingInDegree', 0.0) parameters.put('mapProjection', "WGS84(DD)") terrain = GPF.createProduct('Terrain-Correction', parameters, flood)

Hope this helps!

1 Like

Thank you for response. In my case path is different:

~/snap/bin/gpt -h Terrain-Correction

As I can see, within the output, in “Graph XML Format” section I have:

<demResamplingMethod>string</demResamplingMethod>

which still might be a source of confusion, since you need to know exact spelling of available options (in your example you set BILINEAR_INTERPOLATION). Okay, in previous lines of the output there’s information: Default value is ‘BILINEAR_INTERPOLATION’ but what are the other available options?


Coming back to my question:

I used mentioned tool from command line to retrieve information about ‘Resample’ operator:

~/snap/bin/gpt -h Resample

ResampleOp.txt (2,9 KB)

and I created a stub functions (I assume that resolution is in meters?)

    # I will use Sentinel and SMOS data
    # I will use two functions: getCoarseResProd (to 25 km resolution) or getBetterResProd (to 1 km resolution)

    def getCoarseResProd(file1, destinationPath):
        return getResampled(file1, destinationPath, resolution=25000)

    def getBetterResProd(file1, destinationPath):
        return getResampled(file1, destinationPath, resolution=1000)

    def getResampled(file1, destinationPath, resolution=25000):
        # TODO: this should be tested!!!
        # More info: http://forum.step.esa.int/t/aggregation-and-interpolation-of-sentinel-products-should-i-use-snappy-or-gdal-tools/2522/3
        import snappy
        from snappy import GPF
        from snappy import ProductIO

        product = snappy.ProductIO.readProduct(file1)

        HashMap = jpy.get_type('java.util.HashMap')
        GPF.getDefaultInstance().getOperatorSpiRegistry().loadOperatorSpis()
        parameters = HashMap()
        parameters.put('sourceProduct', product)
        parameters.put('upsampling', "Bilinear")
        parameters.put('downsampling', "Mean")
        # As I checked in SNAP desktop, 'targetResolution' option is sometimes not available
        # and I need to use targetHeight and targetWidth instead
        parameters.put('targetResolution', resolution)
        terrain = GPF.createProduct('Resample', parameters, destinationPath)

        product.dispose()

        return destinationPath

However, I have doubts regarding targetResolution option.
As I check in SNAP desktop I am not always able to use it and I need to provide ‘targetHeight’ (together with targetWidth) option value instead.
However, I don’t understand how to use it

In SNAP desktop, I changed value to much higher (as I understand higher cells should means coarser resolution? In other words, what should I expect if I provide higher value for this option?), but when I saved product… I received 660 GB (Gigabytes!) files instead of previous 250 MB (Megabytes). Both - input and output were saved in .dim (Beam) formats.

@fabricebrito, @marpet Could you please advice?

That’s true the other options should be available too.
For this we either need to extend the parameter annotation with the valueSet setting or, probably even better, change the resamplings methods into an enumeration.
Currently possible values are:
NEAREST_NEIGHBOUR, BILINEAR_INTERPOLATION, CUBIC_CONVOLUTION, BISINC_5_POINT_INTERPOLATION, BISINC_11_POINT_INTERPOLATION, BISINC_21_POINT_INTERPOLATION, BICUBIC_INTERPOLATION

I’ve created two tickets for it: SITBX-396 and SNAP-586

The targetResolution can only be used if the source is reprojected to a map. that’s the reason why you can’t always use this option. The height and width shall give you the option to ensure a certain size of the target.
So I can imagine that it is not easy to find the right settings for data which is not projected to a map. Probably we have to enhance the Resampling Maybe you can reproject first and then do the resampling?

I’ve tried, but without success.
I wrote following function:

def getReprojected(file1, destinationPath, crs=‘EPSG:4326’):
import snappy
from snappy import GPF
from snappy import ProductIO

product = snappy.ProductIO.readProduct(file1)
HashMap = jpy.get_type('java.util.HashMap')
GPF.getDefaultInstance().getOperatorSpiRegistry().loadOperatorSpis()
parameters = HashMap()
parameters.put('source', product)
parameters.put('crs', crs)
parameters.put('resampling', "Nearest")
GPF.createProduct('Reproject', parameters, destinationPath)
product.dispose()

but when I try to execute it, I receive following error message:

ValueError: cannot convert a Python ‘str’ to a Java ‘org.esa.snap.core.datamodel.Product’

Could you please advice?

In the code you have provided I see some problems.
The product should not go into the parameters map. It is provided separately to GPF. The Exception is caused by providing the destinationPath where actually the source product is expected.
To write the result you can either use ProductIO or the Write operator.
So I would change your code to:

  ...
  parameters = HashMap()
  parameters.put('crs', crs)
  parameters.put('resampling', "Nearest")
  result = GPF.createProduct('Reproject', parameters, product)
  # either
  ProductIO.writeProduct(result,  destinationPath, 'BEAM-DIMAP')
  # maybe better regarding performance
  write((result,  destinationPath, 'BEAM-DIMAP'))

def write(product, path, formatName):
  parameters = HashMap()
  parameters.put('file', path)
  parameters.put('formatName', formatName)
  GPF.createProduct('Write', parameters, product)

Thank you very much.

I amended my code as you suggested ( https://github.com/kedziorm/mySNAPscripts/commit/a50dba953d5c3670272293add63687b51041b956 ).

However, now I receive following error message:

org.esa.snap.core.gpf.OperatorException: SAR products should be terrain corrected or ellipsoid corrected

According to https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/product-types-processing-levels/level-1 :

The ellipsoid projection of the GRD products is corrected using the terrain height specified in the product general annotation

so I don’t understand why I need to do correction (why this error does not occur when I use CalibrationOp or SubsetOp?)

Anyway as I can see there are few correction types available in SNAP desktop:

Which one should I choose?

I can’t tell you which one to use. I’m not a SAR expert.
Maybe @lveci or @mengdahl or someone else can help here.

GRD products are not terrain corrected. Due to the nature of SAR acquisitions, in order to get an accurate geocoding you need to account for the geometry of the acquisition. Simply reprojecting will not give you an accurate geocoding. Over scenes where you have a DEM, you should use range Doppler terrain correction. Over ocean or ice you can use average height ellipsoid correction.

1 Like

@lveci and @marpet Thank you for response, I’ve tried to use mentioned operator in getTerrainCorrected function ( https://github.com/kedziorm/mySNAPscripts/commit/e34f40c4d763a337baf2f4d010e596c81d04bb8d )

However, regardless of using it on processed or just downloaded from ESA file, I receive an exception: “java.lang.ArithmeticException: / by zero”

Could you suggest any solution/workaround?

I have also some question regarding this operator (I feel that this might be moved to different topic):

  1. Is it important to apply terrain correction before any other operations (such as CalibrationOp or creating a subset)?

  2. Could you tell me if I should define “demName” parameter as one of the dem names available for download or this will be set automatically?

  3. How “saveDEM” parameter works? My understanding is that if I use it, I will not need to re-download DEM each time when I try to accomplish terrain correction - is it correct?

It seems that after defining additional parameters:

parameters.put(‘demName’, “SRTM 3Sec”)
parameters.put(‘externalDEMApplyEGM’, True)
parameters.put(‘pixelSpacingInMeter’, “0.0”)
parameters.put(‘nodataValueAtSea’, True)

it works properly now

Hello @marpet!

Would you have a code snippet for using the Binning operator in snappy? We’ve ran into some issues for the aggregator part.

Cheers!

There is no example yet. Can you tell me a bit more about the problem you have?

I’m new to Sentinel and SNAP. I’ve looked through my running SNAP and can’t find an option for range Doppler terrain correction. I’m trying to use GDAL 2.4 to read the s1 grd geoTiff image (as downloaded from the hub) and the call to GetProjectionRef() returns an empty string (`’) - gdalinfo does the same thing so I’m pretty sure my code isn’t broken. Because of this I can’t do a ground to image and image to ground mapping. Suggestions very much appreciated.

You can find Range Doppler Terrain Correction in the menu at
Radar / Geometric / Terrain Correction / Range-Doppler Terrain Correction

In the upper right corner of SNAP you can find a little search box. It can be helpful when you are looking for things.