Snappy subsetting problems for S3A data

mariesmith · May 17, 2017, 8:01am

Hi there,

I have been processing data in python with snappy since January 2017 (as part of the S3VT we
receive data for our region from the ftp rolling archive.) However, I have been getting an error for files
downloaded after about 5 April 2017. The last file that I am able to subset successfully within
python is:
S3A_OL_2_WRR____20170405T072830_20170405T081238_20170405T100300_2648_016_149______MAR_O_NR_002.SEN3

However, for files after that date, like this one from the 8th of April:
S3A_OL_2_WRR____20170408T075036_20170408T083446_20170408T102644_2650_016_192______MAR_O_NR_002.SEN3

I have been receiving the following error:

Traceback (most recent call last):
File “”, line 11, in
File “E:\S3A\code\subset_data_test.py”, line 24, in subset_beng
sp = op.getTargetProduct()
RuntimeError: java.lang.IllegalArgumentException: Points of LinearRing do not form a closed linestring

I have a similar error when I try to use Graph Builder for the 8th (Error: Points of LinearRing do not form a closed
linestring), while it still works for dates before that. I include my python snappy subsetting code, which I have
previously applied successfully in python, below. As far as I can see my polygon should be fine?

import snappy
from snappy import ProductIO
import numpy as np
import sys

def subset_beng(fname,filename,options):
SubsetOp = snappy.jpy.get_type(‘org.esa.snap.core.gpf.common.SubsetOp’)
WKTReader = snappy.jpy.get_type(‘com.vividsolutions.jts.io.WKTReader’)

wkt ='POLYGON ((16 -30, 20 -30, 20 -35, 16 -35, 16 -30))'
geom = WKTReader().read(wkt)

p = ProductIO.readProduct(fname+'/xfdumanifest.xml') 
op = SubsetOp()
op.setSourceProduct(p)
print("Getting subset of "+filename[0:24])
op.setGeoRegion(geom)
sp = op.getTargetProduct()
if sp.getSceneRasterWidth()>0:
    print("Writing subset of "+filename[0:24])
    ProductIO.writeProduct(sp, options.subdir+filename+"_beng_subset.nc", "NETCDF4-BEAM")
    print("Done with subsetting")
else:
    print("Area not in orbit. Carry on to next file")

I am still able to subset my region of interest in SNAP manually when exporting to NetCDF4-BEAM
format, and then continue processing in python. However, I need to automate this step for future operational
processing chains.

Has something changed in the file structures? Please help

Kind regards,
Marie

marpet · May 17, 2017, 1:13pm

Hi Marie,

I’ve no access to this archive. Can you upload both products to our ftp?
I’ll send you a separate message with the coordinates.

Have you asked Eumetsat too, if they made any changes?

mariesmith · May 17, 2017, 1:17pm

Hi Marco,

No I have not asked Eumetsat.

I will upload the two files asap

Cheers,
Marie

mariesmith · May 17, 2017, 1:38pm

Files have been uploaded

marpet · May 23, 2017, 10:02am

Hello Marie,

If found the cause for this problem. OLCI RR Orbits can contain the north pole. And this wasn’t properly handled for all cases in our code (SNAP-759). This is fixed and will be released with the next update. Not sure when will be published, maybe in two week or so.

As a temporary workaround, you can do a pixel based subset first and afterwards a geo-region subset.
The pixel based subset might look like this:

from snappy import Rectangle
op1 = SubsetOp()
op1.setSourceProduct§
size = p.getSceneRasterSize()
op1.setRegion(Rectangle(0, 500, size.width, size,height))
pixelSubset = op1.getTargetProduct()

and in your second op you use pixelSubset instead of p

op.setSourceProduct(pixelSubset)

mariesmith · May 23, 2017, 10:29am

Thanks Marco, I will try this out. I did also write a temporary manual subset code in python (although I’m not really happy with it), so I should be able to cope until the next update release. Thank you for finding the problem!

Shaila · April 5, 2022, 8:52am

Hello, I have the same error for SY_2_AOD___ products during a subset with snappy:

RuntimeError: java.lang.IllegalArgumentException: Points of LinearRing do not form a closed linestring

I’m using SNAP 8. I suppose this bug is not fixed yet and I have to subset through pixel or I’m doing something wrong?

here my code:
parameters_subset = snappy.HashMap()

parameters_subset.put('sourceBands', all_band_string[:-1])

parameters_subset.put('copyMetadata','true')

parameters_subset.put('geoRegion','POLYGON(' + AOI_coordinates + ')')

subset =  snappy.GPF.createProduct('Subset', parameters_subset, file_nc)

gnwiii · April 5, 2022, 11:26am

This usually means the last point in the POLYGON is not identical to the first point.

Shaila · April 5, 2022, 12:37pm

The same AOI works with other products, however the output it’s:
POLYGON((12.171435706074094 42.99876354664382,12.171435706074094 43.281325325696436,12.573037857457106 43.281325325696436,12.573037857457106 42.99876354664382,12.171435706074094 42.99876354664382))

I don’t think the problem is the AOI?

gnwiii · April 5, 2022, 1:32pm

So no “north pole” in the polygon. Are any of the points outside the product footprint?

Shaila · April 5, 2022, 4:02pm

nope, no point outside the product. I made this request in the query: footprint:“Intersects(POLYGON((12.171435706074094 42.99876354664382,12.171435706074094 43.281325325696436,12.573037857457106 43.281325325696436,12.573037857457106 42.99876354664382,12.171435706074094 42.99876354664382)))” and i also checked manually on SNAP

gnwiii · April 5, 2022, 6:58pm

Having ruled out the usual suspects, you need to provide enough detail so others can reproduce the problem. You mention that the AOI works for other products – are any of the other products super-pixel data? Did you convert the SY_2_AOD__ products to plain old rasters?

Shaila · April 6, 2022, 10:26am

The main pourpose of my code is to translate Sentinel 3 products in BigTIFF- GeoTIFF data and in order to do so I use snappy to create a subset and then i reproject the result into CRS EPSG:4326 and write the product as GeoTIFF-BigTIFF. This procedure works for ‘SY_2_SYN___’,‘SY_2_VG1___’,‘SY_2_VGP___’ products. Here the part of my code that make the subset and reprojection with snappy:
#SUBSET

    parameters_subset = snappy.HashMap()

    parameters_subset.put('sourceBands', all_band_string[:-1])

    parameters_subset.put('copyMetadata','true')

    parameters_subset.put('geoRegion','POLYGON(' + AOI_coordinates + ')')

    subset =  snappy.GPF.createProduct('Subset', parameters_subset, file_nc)

   

    #RERPOJECT

    parameters_reprojection = snappy.HashMap()

    parameters_reprojection.put('addDeltaBands','false')

    parameters_reprojection.put('crs','EPSG:4326')

    parameters_reprojection.put('includeTiePointGrids','true')

    parameters_reprojection.put('noDataValue','0')

    parameters_reprojection.put('resampling','Nearest')

    reprojection_other_bands = snappy.GPF.createProduct('Reproject', parameters_reprojection, subset)

    output_tmp_allbands = output_tmp + 'file_nc'

    snappy.ProductIO.writeProduct(reprojection_other_bands, output_tmp_allbands, 'GeoTIFF-BigTIFF', loading_bar)

I’m sorry I’m really new in this subject! Have I answered all the questions you made?

gnwiii · April 6, 2022, 12:22pm

You should mention your platform, SNAP Version with update status, and the name of the SY_2_AOD__ file so someone can try to reproduce the issue. Sentinel-3 Spatial Resolution singles out the AOD product:

The SYN AOD products are provided on a wider resolution, created from OLCI image grid. Each super pixel is representing 15 by 15 300m OLCI pixels. Resolution is then 4.5 km.

It is possible that SNAP doesn’t handle this format, or there could be a problem with the particular file. Can you view the SY_2_AOD___ products in the SNAP GUI? If so, you can try adding your polygon to see if it is being shown properly.

Shaila · April 7, 2022, 8:29am

SO:

Edition	Windows 11 Pro
Version	21H2
Installed on	‎08/‎03/‎2022
OS build	22000.556
Experience	Windows Feature Experience Pack 1000.22000.556.0

SNAP version
SNAP 8.0.9

File Used:
https://scihub.copernicus.eu/dhus/odata/v1/Products(‘d3b36308-d55b-4a8b-b7f7-a766d7c8a912’)/$value

Error occurred during the subset creation through snappy ::
java.lang.IllegalArgumentException: Points of LinearRing do not form a closed linestring

The subset creation through SNAP GUI works fine in the same AOI i used in my code, it’s correctly geolocated and conserves fine all the metadata and bands

marpet · April 8, 2022, 7:36am

Do your AOI_coordinates contain parantheses?
Otherwise, there are missing two.

Try:
parameters_subset.put('geoRegion','POLYGON((' + AOI_coordinates + '))')

Shaila · April 8, 2022, 7:49am

My AOI_coordinates has parenthesis, i used this variable also in my query to copernicus scihub. Here an example of query:

(footprint:“Intersects(POLYGON((12.171435706074094 42.99876354664382,12.171435706074094 43.281325325696436,12.573037857457106 43.281325325696436,12.573037857457106 42.99876354664382,12.171435706074094 42.99876354664382)))”) AND ( beginPosition:[2022-02-01T00:00:00.000Z TO 2022-03-04T23:59:59.999Z] AND endPosition:[2022-02-01T00:00:00.000Z TO 2022-03-04T23:59:59.999Z] ) AND ( (platformname:Sentinel-3 AND producttype:SY_2_AOD___) )

in this case the value of AOI_coordinates variable is (12.171435706074094 42.99876354664382,12.171435706074094 43.281325325696436,12.573037857457106 43.281325325696436,12.573037857457106 42.99876354664382,12.171435706074094 42.99876354664382) so i confirm the variable has parenthesis