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:
from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt
from datetime import date
from snappy import GPF,HashMap,jpy
from snappy import (ProgressMonitor, VectorDataNode,
WKTReader, ProductIO, PlainFeatureFactory,
#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
#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 = targetBand1
parameters = HashMap()
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')
#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:
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)
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
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
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.
RuntimeError: java.lang.IllegalStateException: raster data not loaded
While as I pointed before if I do:
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.
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).
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?
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.