Using SubsetOp.getTargetProduct() to clip S3 SLSTR LST data

When using SubsetOp.getTargetProduct() to clip S3 SLSTR LST data an “RuntimeError: org.esa.snap.core.gpf.OperatorException”

Here is a code:

import snappy
from snappy import ProductIO

sen3_in_filename = “H:/sentinel-3/kotelny_island/S3A_SL_2_LST____20180829T031813_20180829T045912_20180830T094318_6059_035_118______LN2_O_NT_003.SEN3”

sen3_out_filename = “H:/sentinel-3/kotelny_island/test.tif”

product = ProductIO.readProduct(sen3_in_filename)

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

lat_min = 73
lat_max = 77
lon_min = 133
lon_max = 152

geom = WKTReader().read(“POLYGON(({} {}, {} {}, {} {}, {} {}, {} {}))”.format(lon_min, lat_min, lon_min, lat_max, lon_max, lat_max, lon_max, lat_min, lon_min, lat_min))

op = SubsetOp()
op.setSourceProduct(product)
op.setGeoRegion(geom)

sub_product = op.getTargetProduct()
ProductIO.writeProduct(sub_product, sen3_out_filename, ‘GeoTIFF’)

then error appears:

sub_product = op.getTargetProduct()
RuntimeError: org.esa.snap.core.gpf.OperatorException

When I use similar code for SMOS L1/L2 data (for example SM_REPB_MIR_SCLF1C_20110801T151944_20110801T161257_620_183_1.DBL). Everything works fine. New raster is created, clipped and saved.

I have read, that maybe it is because of an old version of snappy where was a problem with S3 data when data were close to poles. I have updated SNAP Desktop application. Then create new snappy package using this command:

“C:\Program Files\snap\bin\snappy-conf.bat” “D:\OSGeo4W64\bin\python.exe” “D:\OSGeo4W64\apps\Python27\lib\site-packages”

Just in case I deleted .snap in user directory and snappy folder in python lib directory (in my case “D:\OSGeo4W64\apps\Python27\lib\site-packages”).

When I use SubsetOp without op.setGeoRegion(geom), new raster is created.

Maybe I have still some old version of snappy? When I type ‘pip list’ there is a ‘snappy (2.6.1)’. But I did not use pip to install snappy.

Is there some idea where could be a problem? Or what is the proper way of update snappy?

Zdenek

Well, finally I have successfully updated snappy to version 6.0.4 (as written after ‘pip list’ command). But a problem is still the same…

It looks like the same problem as here - Work with data S3A_SL_2_LST in snap. If I set values of longitude min and max as from image boundary new subset image is created and saved.

Yes, I think it’s the same bug. It will be fixed in the next release.

Do you have an estimate on when the new version will be released. A rough estimate would also be useful?

Hi.
I’m not sure I understood your answer, did you succeed in cropping your LST data ?
Could you provide a code snippet if you indeed did it ?
Thanks a lot

Hi,

no :slight_smile: I did not start to use S3 SLSTR LST data at all. The last what I had tried was cropping file using pixels coordinates. As I have been investigating Kotelny Island, which is quite on the North, I just cropped it by using this rectangle - Rectangle(0, 0, size.width, 2000), where 2000 means, that I have read 2000 rows from the beginning. And it was enough to me. The workaround is about to use the same width as in the original file and get just that rows which contain an area of my interest.

sen3_in_filename = “original_file.SEN3”
sen3_out_filename = “clipped_file.dim”

product = ProductIO.readProduct(sen3_in_filename)

SubsetOp = snappy.jpy.get_type(‘org.esa.snap.core.gpf.common.SubsetOp’)

op1 = SubsetOp()
op1.setSourceProduct(product)
size = product.getSceneRasterSize()
op1.setRegion(Rectangle(0, 0, size.width, 2000))
pixelSubset = op1.getTargetProduct()

# now you can save it
ProductIO.writeProduct(pixelSubset, sen3_out_filename, ‘BEAM-DIMAP’)

Next idea was to use Binning and some aggregation function to get another product, already binned, and then cropped it again with SubsetOp, but now with polygon of lat-lon coordinates (using WKTReader for example).

You can try this for binning:

aggregator_average_config = snappy.jpy.get_type(‘org.esa.snap.binning.aggregators.AggregatorAverage$Config’)

agg_avg_lst = aggregator_average_config(‘LST’)

HashMap = snappy.jpy.get_type(‘java.util.HashMap’)
ArrayList = snappy.jpy.get_type(‘java.util.ArrayList’)

parameters = HashMap()

aggregators = snappy.jpy.array(‘org.esa.snap.binning.aggregators.AggregatorAverage$Config’, 1)
aggregators[0] = agg_avg_lst

parameters.put(‘outputFile’, ‘out.dim’)
parameters.put(‘numRows’, 13360)
parameters.put(‘aggregators’, aggregators)
# pixelSubset is Product created before using SubsetOp
binned = snappy.GPF.createProduct(‘Binning’, parameters, pixelSubset)

and then you can try SubsetOp with lat-lon. But as I mentioned before, I did not try it. Check it and write if it works :slight_smile:

Hey @CzendaZdenda

Thanks a lot for your reply, I did not check if it worked because the problem has been solved and it seems you can now directly use

op = SubsetOp()
op.setSourceProduct(products)
op.setGeoRegion(geometry)
sub_product = op.getTargetProduct()

where geometry is your ROI

Thanks

So I’ve begun to use the subsetOp function now that it works. Unfortunately I don’t seem to have coherent results. the data I’m workin with is the result of the following query :

products = api.query(footprint, date=(‘20181201’, ‘20181202’), producttype=‘SL_2_LST___’, timeliness=‘Non Time Critical’)

And the map I’m using is a rectangle coarsely centered around France (here is a geojson corresponding to the ROI):

{
“type”: “FeatureCollection”,
“features”: [
{
“type”: “Feature”,
“properties”: {},
“geometry”: {
“type”: “Polygon”,
“coordinates”: [[[ -4.8,42.4],[11.5,42.4],[11.5,50.8],[-4.8,50.8],[-4.8,42.4]]]}}]
}

I’m then cropping the whole data to my ROI by doing :

SubsetOp = jpy.get_type(‘org.esa.snap.core.gpf.common.SubsetOp’)
WKTReader = jpy.get_type(‘com.vividsolutions.jts.io.WKTReader’)
footprint = geojson_to_wkt(read_geojson(‘map_test.geojson’))
geometry = WKTReader().read(footprint)
op = SubsetOp()
op.setSourceProduct(products)
op.setGeoRegion(geometry)
sub_product = op.getTargetProduct()

And the result I get is weird, I have nothing but a small region, here is the output area on snap :

@marpet do you happen to have an idea of what’s going on ?

Thanks

Hi, I did not have any follow-up on this issue, I tried reinstalling SNAP and I have the toolboxes up to date.

I re-downloaded data and tried to crop it on my region of interest. What I noticed is the following behaviour :

  • when the ROI is strictly inside the Sentinel 3 data, the crop operation seems to work fine.
  • however if the sentinel 3 does not strictly include the ROI, the behaviour is weird (as depicted above : the ROI corresponds to France and the sentinel 3 data does not cover the whole ROI for this entry).

What I’m intersted in is the temperature over a specific region, and my aim was to download data for each day, then crop to this region each entry of the day and work with this cropped data.

Because of this issue I don’t really know how to handle those data, does anyone have hints for what I should do ?

Thanks again,
Sylvain

How about checking if each corner of the ROI is inside the image. if any corner is out side of the image just discard that image.

Thanks @aksshre, I began by doing this and now I just calculate the intersection of the ROI and the satellite image and make the crop over this intersection, which seems to work nice, so thanks for the hint !

Anyhow I don’t know what is the expected behavior when trying to crop to a ROI which isn’t strictly inside, but I think the actual way it works is a bit buggy

Hi Sylvain,

I investagted your problem earlier but forgot to reply. Sorry.

It seem this is caused by the Geo-Coding we create for the data.
To work around this issue you have three options.
One is what you are already doing. Another would be to shrink the data first e.g. europe (excluding the poles) and the last one would be to rename or remove the geodetic_in.nc file in the *.SEN3 directory.

The intendet behaviour is that your initial approach works.
I’ve added your report to our issue tracker.

Thanks a lot it indeed works when I remove the geodetic_in.nc !

@marpet I have a question about using snappy to export masks(e.g snow mask or cloud mask) from S3 L1 RBT data. It can be done??