Handling sentinel data in python ( working with snappy)

Thank you so much for the links, especially the snappy_flh.py was really helpful, as it shows how to process a sentinel product line by line , and many other things in few lines of code. What I wanted to try also was to display a full image for example with matplotlib but this is probably impossible because of the size, is that the point? anyway your suggestion to consider only subsets of a full image can be the right way. The second approach you described is really interesting but I think it could be useful in a more advanced step of the project I am involved in.
Another question came to my mind while looking at the snappy_flh script: is that possible to create exactly a sentinel-2 product file(I mean with the same structure of an original product, with folders and .xml file) in the same way as it can be done to create for example a geotiff product? The java API is so huge that I do not know where to search for what I need.
Thanks so much so far, for everything

Regarding the display of full images with matplotlib I don’t what limitations it has.
Maybe someone else has a suggestion.

Currently it is not possible to write S2 products in the same format (SAFE) as they are delivered. But we have already thought about.

Thanks for the answer, anyway that is not a problem at the moment, it was basically something I was wondering.
I have another question that is related to the python examples of snappy. Lookind deeply into the snappy_bmaths.py and snappy_ndvi.py I have noticed two different ways to create a product. The first is based on band descriptors and it is created with the GPF module, the second is simply created with the methods of the Product and ProductIO classes. I am trying to understand the reason why two different approaches have been used and one explanation could be that in the second case the creation of specific flags needs a particular processing that can not be implemented just with band maths operations . And probably this could be the same also if we want to create specific masks for a target product. Could you help me understand this aspect?

Both use different approaches. snappy_bmaths.py configures an operator (‘BandMaths’)of GPF. This operator as some specific requirements on its parameters. So here it just uses an existing implementation. You can also invoke such operators from the command line.
snappy_ndvi.py implements the ndvi itself.
What’s best is really depending on your use-case. But if you want to create masks it is probably better to do it more like in the ndvi example.
In the snappy_geo_roi.py example the is an example for a mask based on a shapefile. For shapefiles masks are created automatically when added to the product.

I went to see all the remaining python examples and they were actually really interesting. I have various questions about them but now I will provide the most urgent.
I tried to use the snappy_subset.py script, and as I wanted to avoid problems with memory storage I have tested it with a sample MERIS product that had been provided for BEAM tool. What I did was first to open the MERIS product with SNAP, creating a geometry vector and extract the WKT from the geometry. Then I imported it into the snappy_subset.py and I have to say that it was successful. I tried to do the same with a S-2 product but in some case it is not possible ( when the product has rasters of different sizes) . It seems quite logical, so it means that I could get a subset only for products that have all the rasters with the same size? Another error that I get when trying to obtain a subset from a S-2 product is : RuntimeError: java.lang.ArrayIndexOutOfBoundsException: 0
This comes out for the line where I have the readProduct method and it seems to be related to the array indexes but it is not clear for me.
How could I solve these problems?

Regarding the ArrayIndexOutOfBoundsException I can’t help you. At least I would need some more information. Is more printed as just the exception. Some stack trace below.
You might be right that it might be related to the multi size issue.

Subsets can currently only be done on product where all bands have the same size. In the future this should also work for multi-size products.
Before you Subset you should use the resample operator, to make the bands equal in size.

HashMap = snappy.jpy.get_type('java.util.HashMap')    
snappy.GPF.getDefaultInstance().getOperatorSpiRegistry().loadOperatorSpis()
source = snappy.ProductIO.readProduct('G:/EOData/Meris/RR/MER_RR__1PNBCG20050709_101121_000001802038_00466_17554_0001.N1')

parameters = HashMap()
parameters.put('paramName', False)
result = snappy.GPF.createProduct('Resample', parameters, source)

You can get the available parameters either by calling gpt Resample -h on the command line or by open the Resample Operator in the Desktop App and the go to the menu File/Display Parameters.

Really thanks for the quick replay. It is really kind from you and I am appreciating it so much.
Actually I did not consider the resampling and in my case it is surely what I need to do.
For all the MERIS sample datasets I had no problem while extracting a subset, and none of them has the multi size issue.

I also noticed another element, in this case it is not an issue, but just about the product structure.
When running the snappy_ndvi.py the product that is returned has the FLAG codings node that is really useful to investigate invidually elements like coastline and water. But the Flag codings node and the correspondant masks are not included in the standard S-2 product if I am not wrong. So in theory it would be possible to do it by ourselves?
In my case I think I could get the flags, but what about creating and setting the corresponding masks?
I have found the mask constructor but I do not understand how to put the data(boolean) into it. I wonder if it could be done with a GPF operator as it has been done to create a band.
Thanks in advance again

How to create the flags you see in the snappy_ndvi.py example.
To create corresponding masks I can give you this example:

Color = jpy.get_type('java.awt.Color')
transparency = 0.3
maskColor = Color(200, 100, 100) # Color.YELLOW
targetProduct.addMask('maskName', 'flagcoding.flagName', 'description', maskColor, transparency)

Yes, infact I managed to create flags by myself and I also found the addMask method to create and add a new mask before reading this last reply. But I could not find what should I insert as ‘description’, because the band math requested is included in the 2nd parameter, and so I do not understand what that string should represent. So the product I created actualIy had a new mask but I could not display the corresponding image. I searched for examples or similar questions but I could not see anything.
All these questions come because of the project I am involved in. I can tell you shortly about it, as you are helping me so much…

Hi,

I have tried to implement the resample operator according to the lines of code you suggested, but I got an error that I am going to report here, I also provide the whole script so that you can look at it and maybe it will make clear where is the problem . The script should resample a S-2 product and then extract a subset according to a certain polygon that I obtained using the desktop version of SNAP:

import sys

import snappy
from snappy import ProductIO

SubsetOp = snappy.jpy.get_type(‘org.esa.snap.core.gpf.common.SubsetOp’)
WKTReader = snappy.jpy.get_type(‘com.vividsolutions.jts.io.WKTReader’)

HashMap = snappy.jpy.get_type(‘java.util.HashMap’)
snappy.GPF.getDefaultInstance().getOperatorSpiRegistry().loadOperatorSpis()

parameters = HashMap()
parameters.put(‘paramName’, False)

wkt =‘POLYGON ((11.911873008057473 42.94607650372235, 12.261165250682525 42.93664951680847,
12.253737279745241 42.795528178921806, 11.905237724534322 42.804909016039986, 11.911873008057473 42.94607650372235))’

geom = WKTReader().read(wkt)

print(“Reading…”)
product = ProductIO.readProduct(‘C:\Program Files\snap\S2A_OPER_PRD_MSIL1C_PDMC_20160520T212837_R037_V20160520T112904_20160520T112904.SAFE\S2A_OPER_MTD_SAFL1C_PDMC_20160520T212837_R037_V20160520T112904_20160520T112904.xml’)
result = snappy.GPF.createProduct(‘Resample’, parameters, product)
op = SubsetOp()
op.setSourceProduct(result)
op.setGeoRegion(geom)

sub_product = op.getTargetProduct()

print sub_product.getSceneRasterWidth()

print(“Writing…”)
ProductIO.writeProduct(sub_product, “snappy_subset_output_2.dim”, “BEAM-DIMAP”)

print(“Done.”)

And this is the error report I got :

Traceback (most recent call last):
File “”, line 1, in
File “C:\Anaconda\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py”, line 685, in runfile
execfile(filename, namespace)
File “C:\Anaconda\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py”, line 71, in execfile
exec(compile(scripttext, filename, ‘exec’), glob, loc)
File “C:/Users/antonio/Documents/Python Scripts/Snap/snappy_subset_1_resamp.py”, line 29, in
result = snappy.GPF.createProduct(‘Resample’, parameters, product)
RuntimeError: org.esa.snap.core.gpf.OperatorException: Resample: The Product ‘S2A_OPER_MTD_SAFL1C_PDMC_20160520T212837_R037_V20160520T112904_20160520T112904’ already contains a band with the name ‘Intensity_null’.

What could be the reason? should I have given some parameters setting for the Resample operator? I thought that using

parameters.put(‘paramName’, False)

meant to set the default values of the resample operator, but probably that is wrong.

Thinking again about what you told me about the resampling operation I tried to look at the new SNAP 4.0 after installing it on another PC and I discovered that a new resample operator has been added besides the one that was in the RADAR menu and the resampling options that appear when opening a Sentinel product. So probably when you suggested how to use the resampling operator you were referring to the newest version because I am working with SNAP 2.0 and the desktop version has a differet resample menu. Is that correct?

I also have another question that derives from the snappy_write_image.py. I tried on my own to save a gray-scaled image from a band and the RGB image obtained from 3 bands. While in the example the image had been saved as a RenderedImage object it seems to me that I had to use the BufferedImage class for my purpose. I am going to post the report of my attempt for the RGB image, but the result was saving a png white image instead of the one I wanted . I tried to do it from a Python console, not with a script.

from snappy import jpy
from snappy import ProductIO
from snappy import ProductUtils
from snappy import ProgressMonitor
image_info = jpy.get_type(‘org.esa.snap.core.datamodel.ImageInfo’)
data = jpy.get_type(‘org.esa.snap.core.datamodel.RGBChannelDef’)
imageio = jpy.get_type(javax.imageio.ImageIO)
file = jpy.get_type(‘java.io.File’)

p = ProductIO.readProduct(‘C:\Program Files\snap\MER_RR__1PNPDK20030813_175754_000026132019_00027_07596_4557.N1’)

BLU = p.getBand(‘radiance_1’)
GREEN = p.getBand(‘radiance_5’)
RED = p.getBand(‘radiance_13’)
info = data([‘radiance_13’, ‘radiance_5’, ‘radiance_1’])
INFO = image_info(info)
image = ProductUtils.createRgbImage([RED, GREEN, BLU], INFO, ProgressMonitor.NULL)
output = file(‘MINE_RGB.png’)
imageio.write(image,“png”,output)

I did not get any error, but the image was totally white, so I guess I missed some color setting but honestly I do not know how to fix it.

As description you can use any string which describes your mask. Or leave it empty.
But sometimes it can be useful for users to know a bit more than simply knowing that it is a e.g. a land mask.
But the description shouldn’t be to long.

Why the new mask in your case is not displayed I can’t say. I would need some more information to do so.

Yes, I was referring to the Resample operator which was introduced in SNAP 3.

Regarding the image creation I will answer in the separate threads you have created.

Really thanks!! I will work more on the mask description then! Talking about the images,I am going to update the thread because I found some partial solution.
Thanks again for the moment!

Hello Marpet and Antonio,

One question, I need apply resample and subset a Sentinel 2 product.

I have one problem in the parameters of function:
result = snappy.GPF.createProduct(‘Resample’, parameters, product)

The error is:
RuntimeError: org.esa.snap.core.gpf.OperatorException: Either referenceBandName or targetResolution or targetWidth together with targetHeight must be set.

You can help me?

Thanks.
Mobil. +52(442)1770038

If you type on the command line gpt Resample -h you will get the help for the operator.

The operator needs to know to which resolution the data shall be resampled, therefore it is necessary to specify at least one of the parameter.

So you can do one of the following:

parameters.put('referenceBandName', 'B1') # B1 = 60m 

or

parameters.put('targetResolution', '15') # 15m

or

parameters.put('targetWidth ', '2500') # fixed size; resolution is computed
parameters.put('targetHeight ', '2500') # fixed size; resolution is computed

or

Dear all,

I have to assemble different slice more than 70 images for my study area.
So, is it possible to use snappy in python to process lots of images for assembling?
I have configured Python to the SNAP-Python (snappy) interface, but I have no idea how to set the parameters in “Slice Assembly” function in snappy to apply.

The script I have write is as follow:

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

_path1 = “L:\Sentinel images\track_69\test\1S1A_IW_SLC_1SSV_20141022T100021_20141022T100051_002941_003570_A0CE.zip”

sentinel1 = ProductIO.readProduct(path1)

_path2 = “L:\Sentinel images\track_69\test\2S1A_IW_SLC_1SSV_20141022T100049_20141022T100119_002941_003570_8A60.zip”

sentinel2 = ProductIO.readProduct(path2)

GPF.getDefaultInstance().getOperatorSpiRegistry().loadOperatorSpis()

_HashMap = jpy.get_type(‘java.util.HashMap’)
_SliceAssembly = jpy.get_type(‘org.esa.s1tbx.sentinel1.gpf.SliceAssemblyOp’)

op = SliceAssembly()
_op.setsourceProducts(sentinel_1)
_op.setsourceProducts(sentinel_2)
op.setselectedPolarisations(‘VV’)

product = op.getTargetProduct()

And it goes wrong with "AttributeError: ‘org.esa.s1tbx.sentinel1.gpf.SliceAssemblyOp’ object has no attribute ‘setsourceProducts’"

Please help me or give me some hints.

Many thanks!

You should call GPF to invoke the operator. After reading the products do:

from snappy import HashMap

parameters = HashMap()
parameters.put('selectedPolarisations', 'hh,vv') # or whatever is necessary here
products = jpy.array('org.esa.snap.core.datamodel.Product', 2)
products[0] = sentinel1 
products[1] = sentinel2 
result = GPF.createProduct('SliceAssembly', parameters, products)
ProductIO.writeProduct(result, 'C:\data\temp\result.nc', 'NetCDF4-CF')
1 Like

This is really interesting. I have never needed such structure so far…but I will take note…for future tutorials :wink:

Hi @marpet,

Sorry for disturbing you again! I have tried the code you mentioned, but the result.dim only included Metadata, Vector Data and Tie-Point Grids. The Bands were disappeared.

I’m new in using snappy, so I don’t really understand how to deal with this problem.

My script now is as below:

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

images\track_69\test\S1A_IW_SLC__1SSV_20141022T100021_20141022T100051_002941_003570_A0CE.zip"
sentinel_1 = ProductIO.readProduct(path_1)

path_2 = “L:\Sentinel images\track_69\test\S1A_IW_SLC__1SSV_20141022T100049_20141022T100119_002941_003570_8A60.zip”
sentinel_2 = ProductIO.readProduct(path_2)

from snappy import HashMap

parameters = HashMap()
parameters.put(‘selectedPolarisations’, ‘hh,vv’) # or whatever is necessary here
products = jpy.array(‘org.esa.snap.core.datamodel.Product’, 2)
products[0] = sentinel_1
products[1] = sentinel_2
result = GPF.createProduct(‘SliceAssembly’, parameters, products)
ProductIO.writeProduct(result, ‘L:\Sentinel images\track_69\result.dim’, ‘BEAM-DIMAP’)

Many thanks!

Hi @marpet,

I have solved this problem!
I changed the Java_max_mem to 12G in snappy.ini, and everything goes fine!
Now, I could do SliceAssemmbly of my images automatically
by python!
Thank you very very much!!

1 Like