How to create masked bands from imported WKT geometry in Snappy

Within the snappy framework, I would like to apply a mask to a sentinel-1 image using a wkt string.

I know this can be done with a Land-Sea-Mask, but the geometry needs to be added to the product first.
Adding the geometry can be done directly using the Import-Vector, but it only allows a shp. from file (please correct me if i’m wrong).
To use a wkt string I have tried the example script in snappy_geo_roi.py
(see below for example to try-need to change the IO file path of course, and the file name is original so people want they can search for it here: https://scihub.copernicus.eu/dhus/#/home)

However, the vector it creates is empty (see attached pic) . So when I try and run the Land-Sea-Mask I get an error:
“RuntimeError: org.esa.snap.core.gpf.OperatorException: expression: Undefined symbol ‘shape’. due to Undefined symbol ‘shape’.”

I know this error also appears if the image is out of bounds, but im confident it is within the extent (I drew it manually in the toolbox):

As a temporary workaround, I dont mind if i have to write the wkt to a shapefile, and then import it via ‘Import-Vector’ operation (export to shapefile is possible from the he Sentinel GUI).

If there is another way to import the wkt, please let me know! thanks in advance!

Test script:

> import sys
> import os

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


> test_folder = r'C:\Users\myName\Documents\sentinel_1'                      
> inputFileName = 'S1B_IW_GRDH_1SDV_20161205T061407_20161205T061432_003257_0058D8_C74B.zip'                          
> wellKnownText = 'POLYGON ((-0.3788970856103323 51.80800809372078, -0.376076940574716 51.80800809372078, 
>    -0.376076940574716 51.80640687188327, -0.3788970856103323 51.80640687188327, 
>    -0.3788970856103323 51.80800809372078)) '
> geom_name = 'shape'

> full_path_input = os.path.join(test_folder,inputFileName)

> geometry = WKTReader().read(wellKnownText)

> product = ProductIO.readProduct(full_path_input)

> wktFeatureType = PlainFeatureFactory.createDefaultFeatureType(DefaultGeographicCRS.WGS84)
> featureBuilder = SimpleFeatureBuilder(wktFeatureType)
> wktFeature = featureBuilder.buildFeature(geom_name)
> wktFeature.setDefaultGeometry(geometry)

> newCollection = ListFeatureCollection(wktFeatureType)
> newCollection.add(wktFeature)

> #productFeatures = FeatureUtils.clipFeatureCollectionToProductBounds(newCollection, product, None, ProgressMonitor.NULL)

> node = VectorDataNode(geom_name, wktFeatureType)
> product.getVectorDataGroup().add(node)

> #vdGroup = product.getVectorDataGroup()
> #for i in range(vdGroup.getNodeCount()):
> #    print('Vector data = ', vdGroup.get(i).getName())
> #
> #maskGroup = product.getMaskGroup()
> #for i in range(maskGroup.getNodeCount()):
> #    print('Mask = ', maskGroup.get(i).getName())

> target_name = 'test_empty_geom.dim'

> """Can manually check the outout file in sentinel toolbox for missing 'wkt coords'"""
> target_name = 'test_empty_geom.dim'
> full_path_vector = os.path.join(test_folder,target_name)
> ProductIO.writeProduct(product, full_path_vector, 'BEAM-DIMAP')

> print('Done with test to add product vector')

1 Like

There is something essential missing in the example.
Adding the features to the VectorDataNode.
I’ve updated it:

Thanks for the super fast response! I’ll try it now.

As an aside, is there a docs for the python methods that I have overlooked?

No there is no special documentation for the python methods.
You have to check the Java documentation.
http://step.esa.int/docs/v5.0/apidoc/engine/

Thanks, the mask is created now thanks! … however it creates virtual bands tmpbands which is doubling my bandcount! I would ideally like a product with just the masked pixels (and preferably noData = numpy.nan)

UPDATE:
I am applying the Land-Sea-Mask with ‘bypass = True’ to just create a polygon mask in the mask group node (see full code below).

I have read Snappy: Mask (Clouds and others)
and readValidMask and isPixelValid methods

from:[quote=“marpet, post:8, topic:1654”]
When you call getValidMaskImage() on the result of product.getMaskGroup().get(‘opaque_clouds_10m’) then you retrieve the valid mask for the cirrus mask.
[/quote]

I then try and fetch the mask using:
land_sea_mask = land_sea_mask_product.getMaskGroup().get(str(geom_name)) land_sea_mask.getValidMaskImage()

I get the error:
AttributeError: 'org.esa.snap.core.datamodel.ProductNode' object has no attribute 'getValidMaskImage'
i think because getValidMaskImage() works on org.esa.snap.core.datamodel

@marpet Can you help tell me what I am missing to apply this mask to bands in a for loop?
Thanks in advance!

I currently have:

geom_name = 'shape'
product = ProductIO.readProduct('my_s1_dim_here')  
polygons_wkt = "POLYGON ((-0.379473496337121 51.806829421554639...)"



#add the wkt geometry to the product using a seperate method which returns a product with added vector
product_with_vector = self.add_wkt_vector_geom(product , geom_name, polygons_wkt)

#create the mask in the product        
params = HashMap()
#        parameters.put('landMask', False)
params.put('useSRTM', SRTM)
params.put('geometry', geom_name)
params.put('invertGeometry', invert)
params.put('byPass', bypass) 

#create a mask inside a product using Land-Sea-Mask (bypass = True)
land_sea_mask_product = GPF.createProduct('Land-Sea-Mask', params, product_with_vector) 

h = land_sea_mask_product.getRasterHeight()
w = land_sea_mask_product.getRasterWidth() 
bands = list(land_sea_mask_product.getBandNames())            
k = int(len(land_sea_mask_product.getBands()))        

#retrieve the mask
land_sea_mask = land_sea_mask_product.getMaskGroup().get(str(geom_name))        
valid_mask = land_sea_mask.getValidMaskImage() 

#create array for mask (this probably doesnt need ot be remade with each loop!)
mask_array  = np.zeros(w, dtype=np.uint8)         
       
#make new product and set writer
mask_product = Product('masked', 'masked', w, h)
ProductUtils.copyGeoCoding(product , mask_product) 
writer = ProductIO.getProductWriter('BEAM-DIMAP')
mask_product.setProductWriter(writer)
mask_product.writeHeader(String('mask_from_wkt.dim'))   

#apply to all bands
for band_index in range(k):
    print('reading band index (k): '+str(band_index+1))
    
    #add a band to be masked
    masked_band = mask_product.addBand(bands[band_index]+'masked', ProductData.TYPE_FLOAT32)
    #set the band noData values
    masked_band.setNoDataValue(np.nan)
    masked_band.setNoDataValueUsed(True)
    
    #array for mask values
    masked_array = np.zeros(shape=(w, h), dtype=np.float32)                    
    #load the band to read
    load_band = land_sea_mask_product.getBandAt(band_index) 
    #read the pixels into an array
    load_band.readPixels(0, 0, w, h, masked_array)
                
    
    #read to the the mask
    valid_mask.readValidMask(0, 0, w, h, mask_array) 
#           #create a numpy mask condition            
    invalid_mask = np.where(mask_array == 0, 1, 0)
    #apply the mask to our masked array
    numpy_mask = np.ma.array(masked_array, mask=invalid_mask, fill_value=np.nan)
    #write the masked band
    masked_band.writePixels(0, 0, w, h, numpy_mask)

print('Masking finished')

return mask_product

OK, my first observation is that you are not working with the most recent SNAP version. The ‘byPass’ parameter doesn’t exist anymore. Please check for updates. In the menu at Tools / Plugins.
Second, a mask itself has no valid mask. Only bands have. The mask can be used as valid mask of a band.
Third, writing the product as you try it will not work. Setting the Product writer to the output product is missing.
See

Beside this, I think you don’t need the Land-Sea-Mask operator at all.
Just after reading the geometry you can use the mask already.
I’ve attached a script. I think it should work, even I haven’t let it run.
Maybe it needs some minor tweaks.

writeNans.py (2.1 KB)

Great, I will give it a shot thanks. Also will get on the updates.

Sorry I haven’t shown it, but I have a separate read-write class. I put an quick input method at the beginning for readability for wider forum users, but forgot to do the same for writing…i’ll update it with your addition.
Iv had a look at your script, and the shortened approach make a lot of sense!

Update:
@marpet, I am trying to run it, and coming into the same problem as previously when reading the mask into the mask array in the line:
geom_name.readPixels(0, 0, w, h, mask_array)

I figured this is a small mistake as the variable from:
geom_mask = src_product.getMaskGroup().get(geom_name)

was unused, and readPixels (to mask array) was being called on the string ‘geom_name’. So I tried changing it to:
geom_name.readPixels(0, 0, w, h, mask_array)
but got the error:
AttributeError: 'org.esa.snap.core.datamodel.ProductNode' object has no attribute 'readPixels'

so i tried using the readValidMask method as in the other forum topics:
geom_mask.readValidMask(0, 0, w, h, mask_array)

But still get the same error:
AttributeError: 'org.esa.snap.core.datamodel.ProductNode' object has no attribute 'readValidMask'

Update:
I understand this is because readValidMask() is applied to bands(not masks) as you stated earlier. In which case, how does one read the mask into an numpy array?

As an aside, what band values does readValidMask() return?

Thanks!

Hi @marpet , thanks for your patience and help on this. I am still trying different approaches as couldn’t get script to read the mask into an array.
I have tried

So i tried:

geom_mask = src_product.getMaskGroup().get(geom_name)
valid_mask_image = geom_mask.getValidMaskImage()

and got the error:
AttributeError: 'org.esa.snap.core.datamodel.ProductNode' object has no attribute 'getValidMaskImage'

Given that i need to call readValidMask() on a band and not a mask, I have then tried your following advice:

I wrote it as (after mask_array, and before the invalid_mask):

        src_band.setValidPixelExpression(str(geom_name) +'== 0')
        src_band.setNoDataValue(np.nan)
        src_band.setNoDataValueUsed(True)
        
        src_band.getValidMaskImage()
    
        #read the mask pixels into an mask array
        src_band.readValidMask(0, 0, w, h, mask_array)

I get the error:

 src_band.readValidMask(0, 0, w, h, mask_array)
RuntimeError: no matching Java method overloads found

Also, in this approach, I cant see how to use the geom_mask, or is it even necessary?
geom_mask = src_product.getMaskGroup().get(geom_name)

Thanks again!

P.S. OK that I also changed the name of the thread to make it more specific to the problem?

Very much appreciated! :slight_smile:

Calling getValidMaskImage() on a mask is not working, because a mask doesn’t have a mask itself.

src_band.readValidMask(0, 0, w, h, mask_array)

This should actually work. It depends, are w and h integer values? Do you create mask_array like this?

mask_array = np.zeros(shape=(w, h), dtype=np.uint8)

What might help you is the debug option.
Have a look into the snappy folder and open the __init__.py.
Change line 54 to

debug = True

Then you should see more information about the error.

Hi Yes, the code is much as you wrote it in the script you uploaded-- ‘writeNans.py’ + recommended edits.
mask array is:
mask_array = np.zeros(shape=(w, h), dtype=np.int32)

I have uploaded it the script again with the above edits iv tried, and also added my missing method ‘add_wkt_vector_geom()’ to make it easier for others to follow later. This should run stand-alone now.

The line of interest are: 101 -117

After adding debug = True, I get the following parameters which look OK
(I was testing here with a stacked file preprocessed to gamma):

('Num features = ', 1L)ipdb> 
('Vector data = ', 'pins')
('Vector data = ', 'ground_control_points')
('Vector data = ', 'shape')
('Mask = ', 'shape')
('defining output bands:', 'Gamma0_VV_mst_05Dec2016')
('defining output bands:', 'Gamma0_VV_slv1_11Dec2016')
('reading band:', 'Gamma0_VV_mst_05Dec2016')

mask_with_wkt.py (5.7 KB)

Try changing from int32 to uint8

@marpet Nearly! I have just checked the file and all the pixels put out to 0.00…

could it be this expression?
src_band.setValidPixelExpression(str(geom_name) +'== 0')

Also I can seem to write the out_product using the usual:
ProductIO.writeProduct(out_product, ‘dir_string’, ‘BEAM-DIMAP’)

RuntimeError: java.lang.IllegalStateException: no product reader for band 'Gamma0_VV_mst_05Dec2016'

When do you get the IllegalStateException? I assume in your very last line.

ProductIO.writeProduct(masked_product, target, 'BEAM-DIMAP') 

Just remove it. It is already written.
Provide the target path as parameter to your method and pass it to product.writeHeader().

In general I would like to suggest, that you consider implementing this as an operator. The implementation will be easier if once understood the concept of operators. In addition it will scale better to bigger scenes.
Have a look here:
https://senbox.atlassian.net/wiki/display/SNAP/How+to+write+a+processor+in+Python

Cool thanks for the info! Unfortunately I am still getting all pixel values turning out as zero after the mask.

This method is to form part of a larger workflow {raw product ->prep(orbit/subset bbox/radio calib etc.) -> stack -> get AOI (subset+mask)-> build model}. I am passing products along the chain so just want to return the masked product.

Thanks!

@Marpet thanks for all your help. I ended up using a similar solution to the ndvi_with_masks.py, specifically using the readValidMask() method.

For those following, the key part:

    #apply the mask as an expression
    src_band.setValidPixelExpression(str(geom_name)+' == 0')
    #read the pixels (to be masked) into the array
    src_band.readPixels(0, 0, w, h, data_array)
    #read the mask into the mask array
    src_band.readValidMask(0, 0, w, h, mask_array)

    #Here set the pixels instead of writing them (so we can return a product instead of writing)
    out_band.ensureRasterData()
    out_band.setPixels(0, 0, w, h, mask.filled(np.

One last question…it is possible to apply multi-part polygons as an expression? Instead of masking the image repeatedly for every polygon to get each of the polygons summary statistics, it would be more ideal to extract the values for parts in a multi-part polygon in one go…

I’m glad that you get working.

Yes, I think also a mulit-polygon should be readable. But you can’t distinguish the different polygons
in the mask. But you can read all the polygons, create different vector data nodes and masks and use them in your iteration. If this is what you’re asking, I’m not sure if I got you right.

1 Like

Hi, i am trying to cut and .safe file into an ROI file because the original .SAFE is so big for processing and i just need a small area but i get this error , thanks if you can help me

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


#if len(sys.argv) != 3:
#    print("usage: %s <inputProduct> <wkt>" % sys.argv[0])
#    sys.exit(1)

#input = sys.argv[1]
#inputFileName = input[input.rfind('/')+1:]
#WKTReader = snappy.jpy.get_type('com.vividsolutions.jts.io.WKTReader')
wellKnownText = 'POLYGON ((-74.911862 4.162501, -74.915142 4.162332, -74.914788 4.156032, -74.912351 4.156068))'
file = ('C:\Users\alejandro vergara\Desktop\S2A_OPER_PRD_MSIL1C_PDMC_20160131T020509_R025_V20160130T153109_20160130T153109.SAFE')

geometry = WKTReader().read(wellKnownText)

product = ProductIO.readProduct(file)

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

newCollection = ListFeatureCollection(wktFeatureType)
newCollection.add(wktFeature)

productFeatures = FeatureUtils.clipFeatureCollectionToProductBounds(newCollection, product, None, ProgressMonitor.NULL)

node = VectorDataNode('shape', productFeatures)
print ('Num features = ', node.getFeatureCollection().size())
product.getVectorDataGroup().add(node)

vdGroup = product.getVectorDataGroup()
for i in range(vdGroup.getNodeCount()):
    print('Vector data = ', vdGroup.get(i).getName())

maskGroup = product.getMaskGroup()
for i in range(maskGroup.getNodeCount()):
    print('Mask = ', maskGroup.get(i).getName())

the error is:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-1-e48baaa4177a> in <module>()
     16 file = ('C:\Users\alejandro vergara\Desktop\S2A_OPER_PRD_MSIL1C_PDMC_20160131T020509_R025_V20160130T153109_20160130T153109.SAFE')
     17 
---> 18 geometry = WKTReader().read(wellKnownText)
     19 
     20 product = ProductIO.readProduct(file)

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

For the WKT polygon you need to repeat the first coordinate as last coordinate, otherwise it is not closed.
Probably this solves already the issue.