Snappy ROI problem

Hi there,
I’m practicing a bit with snappy. My final objective is to automatically extract vegetation indices from specific regions (the ones belonging to viticulture fields).
Right now I am stuck with the sub setting of the data of specific regions.
Here is my code:

import snappy
import sentinelsat
from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt
import subprocess
from datetime import date

from snappy import GPF,HashMap,jpy

from snappy import (ProgressMonitor, VectorDataNode,
                    WKTReader, ProductIO, PlainFeatureFactory,
                    SimpleFeatureBuilder, DefaultGeographicCRS,
                    ListFeatureCollection, FeatureUtils)

#Reading an example geojson file I created selecting an area thanks to https://geojson.io:
footprint = geojson_to_wkt(read_geojson('./map.geojson'))

#Passage found in the snappy example on ROI, not sure if i can just skip it and use the just created variable footprint instead:

geometry = WKTReader().read(footprint)

#Reading the product I downloaded previously  thanks to sentinelsat. The product contains the region i'm interested in

product = 
snappy.ProductIO.readProduct('S2A_MSIL1C_20210225T102021_N0209_R065_T32TNR_20210225T123254.zip')

#Here I just create another product called result with only one band to avoid the bug concerning the different resolutions. (I had troubles when I was trying to subset small regions).

BandDescriptor = jpy.get_type('org.esa.snap.core.gpf.common.BandMathsOp$BandDescriptor')

targetBand1 = BandDescriptor()
targetBand1.name = 'B2'
targetBand1.type = 'float32'
targetBand1.expression = 'B2'


targetBands = jpy.array('org.esa.snap.core.gpf.common.BandMathsOp$BandDescriptor', 1)
targetBands[0] = targetBand1

parameters = HashMap()
parameters.put('targetBands', targetBands)

result = GPF.createProduct('BandMaths', parameters, product)

#From here I tried creating the ROI following the example provided in the documentation

wktFeatureType = PlainFeatureFactory.createDefaultFeatureType(DefaultGeographicCRS.WGS84)
featureBuilder = SimpleFeatureBuilder(wktFeatureType)
wktFeature = featureBuilder.buildFeature('shape')
wktFeature.setDefaultGeometryProperty(geometry)

#Here I get the following error: RuntimeError: no matching Java method overloads found. 
# In the example it is wtkFeature.SetDefaultGeometry instead of 
# wtkFeature.SetDefaultGeometryProperty but if I use the former it throws an error saying 
# 'org.opengis.feature.Feature' object has no attribute 'setDefaultGeometry' 

What am I missing?
Furthermore, concerning the resolution bug, I read that I should resample all the bands at the same resolution but how should I do it?
Finally, before following the example I was trying to get the ROI region using:

parameters = snappy.HashMap() 
parameters.put("copyMetadata", True)
parameters.put("geoRegion", geometry)
subset = snappy.GPF.createProduct("Subset", parameters, result)
snappy.ProductIO.writeProduct(subset, "subset", "GeoTiff")

What’s wrong with this approach? I think I am really confused about the conversion between pixels and geographical coordinates. Any kind of help would be super appreciated =D

You have a process with several steps and several problems, so my response will focus on general approaches to multi-step workflows. It is often best to break the workflow into small steps and verify each one in turn. In some workflows there are steps (such as defining polygonal regions) that only need to be done once and may benefit from using formats with good metadata and application support (NetCDF over GeoTIFF, shapefiles over WKT). When you use small steps and encounter problems it is much easier to provide a self-contained example suitable for posting to the forum.

Starting with geojson to WKT, you should be able to save the WKT and view it in the GUI.

You can also work on the resolutions as a separate step. It is often helpful to s
tart with the GUI and then moving to a snappy script. There have been a number of forum posts dealing with this step.

It may also help to use “try except” in your Python script to get better information when errors occur. When you are processing many images you can encounter problems with corrupt files, running out of memory, etc. that are much easier to handle if you have good error messages or log files.

When referring to documentation it is helpful to provide a link, especially if the documentation is outdated and needs correction.

I don’t know why your wktFeature seems to be of type Feature, where it actually should be of type SimpleFeature.
You can try the following:

# create the SimpleFeature type at the beginning of your script once
SimpleFeature = jpy.get_type('org.opengis.feature.simple.SimpleFeature)

# later cast your wktFeature and set the value
simpleFeature = jpy.cast(wktFeature , SimpleFeature)
simpleFeature.setDefaultGeometry(geometry)

A similar issue happend here:

Thanks a lot marpet! That worked!

Now I’m encountering another problem though which is the one described in this other thread.

It looks like the mask doesn’t apply and data is full of just zeros.

mymask = maskGroup.get(‘shape’)
my_real_mask = jpy.cast(mymask, Mask)
h = product.getSceneRasterHeight()
w = product.getSceneRasterWidth()
data = np.zeros(w * h, np.int32)
data.shape = h,w
my_real_mask.readPixels(0,0,w,h, data)

Any clues? =(
Furthermore, when using subset instead of the mask how can I specify that all the pixels not belonging to the geometry provided will be put to zero? The geometry of the subset has to be a rectangle (Am I wrong?). I would like to get 0 on all the pixels inside the rectangle that don’t belong to the geometry provided.
From this link it looks like you already provided this feature “User should have the option to set the pixels outside of the region to NaN values” unless I’m misunderstanding it.
Thanks a lot for your avaibility, I’m really going crazy here :stuck_out_tongue:

Thanks for your kind answer, will try to be more consistent in the future!
I tried to export the bands in the specified region with the GUI and it worked pretty easily! Now I have to automate the process.

The type of subsetting is not yet implemented in the subset operator yet.
But you can first subset your data to the bounding box and afterwards.
Check these posts:

Why the mask provides no data I don’t understand.
Have you tried this?

my_real_mask.getPixels(0,0,w,h, data)

After your suggestion I tried the

my_real_mask.getPixels(0,0,w,h, data)

but i get the following error:

RuntimeError: java.lang.IllegalStateException: raster data not loaded

While as I pointed before if I do:

my_real_mask.readPixels(0,0,w,h, data)

I get no errors but data is just full of zeros.
I doubled checked that the wkt region is inside the downloaded product trough the GUI.
Concerning the Subset with .shp I will check it tomorrow morning and will update you on how it works.

Just a little note, in your answer the last link “Snappy ROI problem” points to this post. Don’t know if it was intended to be another link or it just popped there by mistake.

Thanks again for your fast answer =D

Glad you are making progress.

You should be able to use the GUI to create “mask” bands from the shapefile that can later be used with band-maths to set pixels outside the region to a “no data” value. I should also mention that if you end goal is areal statistics for the regions, then you might be to use StatisticsOP.

I wonder if some of what you are doing is more in the domain of a GIS package such as QGIS or ArcGIS (both support Python scripts).

Okay, could be that you need to call.

my_real_mask.loadRasterData()

before calling getPixels.

If this is not working
try:

maskImage = my_real_mask.getSourceImage()
maskImage.getData().getPixels(0,0,w,h, data)

Tried both

my_real_mask.loadRasterData()
my_real_mask.readPixels(0,0,w,h, data)

and

maskImage = my_real_mask.getSourceImage()
maskImage.getData().getPixels(0,0,w,h, data)

the code doesn’t throw any error but the mask is still full of zeros…

Concerning the Subset I checked the links you provided me but I couldn’t find a method to obtain what I am looking for:
I can extract the bounding region but I can’t figure how to put at zeros all the elements not belonging to the geometry provided. Am I missing something?

Could you try loadRasterData with getPixels()?

Maybe there is some issue with masks. But I can’t imagine what it is at the moment.

The Sea-land-mask operator does net set to zero but sets the not data value.
I’ve used it in the GUI. I think this is the result you would like to achieve.

Sorry for not being precise, I tried all the possible combinations includind this:

my_real_mask.loadRasterData()
my_real_mask.getPixels(0,0,w,h, data)

but still didn’t work.
Could it be a problem with the snappy version as pointed out by this thread?

Concerning the Sea-land mask it’s exactly what i’d like to achieve, thanks! Is there any example on how I should use it with python instead of the GUI?

Maybe, but I don’t see a reason why it still works in SNAP (Java) but it does not in snappy (Python).
However, I’ve created an issue for this, so that someone can have a look at it.
https://senbox.atlassian.net/browse/SNAP-1413.