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
Aggregation and Interpolation of Sentinel products - should I use snappy or GDAL tools?
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)
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 ‘
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:
but I can read here only that ‘setParameter’ method is inherited from Operator class.
Still no idea how to set parameters for Resampling.
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!
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:
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’):
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)
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.
@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):
Is it important to apply terrain correction before any other operations (such as CalibrationOp or creating a subset)?
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?
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”)
it works properly now
Would you have a code snippet for using the Binning operator in snappy? We’ve ran into some issues for the aggregator part.
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.