I have recently installed SNAP in my Windows8.1 with 8GB RAM and as I said in my first message I am trying to work with SNAP in python using Snappy. I managed to use it properly with the test dataset that is provided with the snappy module in the SNAP directory. I usually work with numpy to handle and process images or multidimensional arrays. But the problem is that the test dataset has really small size if compared with the usual sizes of a Sentinel dataset. And when I try to read for example a Sentinel-2 dataset of small size(1-2 GB) , I manage to read it properly but then if I try to store even only one band in a numpy array I get a MemoryError. So I tried to do the same operation with a smaller subset of the image for a single band and that was successful. So what could I do to work with full images in Python without facing those memory problems? Thanks for the reply!
The question is if you really need the full image data or if you could process the data in frames. Iâve mentioned some sources of python examples in this post. Maybe you can get some inspiration there.
The point is that according to our project we would like to create a set of Python script to show how to work with snappy for a case of study that can be valuable for educational purposes, or as a training example. By " process in frames" do you mean to process a full image after a segmentation in various subsets with smaller sizes? in that case would it possible to use the SNAP API on those subset as if they were original sentinel images? I guess I should âsaveâ those subset in the Sentinel-2 format using Snappy . Is that the idea?
You donât need to do the segmentation upfront.
It depends a bit on your use case how to do it.
Mainly there are two ways of processing the data.
One is by simply using the SNAP API.
You can find an example in snappy_flh.py. This file is also located in the examples folder in the snappy directory.
In this file you can see that the source data is read line by line. But you can change it to any rectangular shape. This is what I meant by frame. In SNAP terminology frame is in most cases named tile.
The other approach is to implement an operator.
A simple example can be found in the ndvi_op.py. A real life example can be found at the S2-RUT tool.
In order to develop such an operator in python three methods need to be implemented.
def initialize(self, context):
def dispose(self, context):
and either
def computeTileStack(self, context, target_tiles, target_rectangle):
or
def computeTile(self, context, band, tile):
The computeTile method is preferred and should be used if results can be computed independently. The computeTileStack should be implemented if the resulting bands are all computed at once, e.g. as output of a neural net.
The computeTile method is also preferable for S2 data because this allows to handle the different resolutions. If the computeTileStack is implemented the source needs to be resampled upfront.
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!