Change in polygon size when reading a WKT with snappy

starting from a Level 1C product opened in Sentinel-2 Toolbox ver 8.0, I choose the B4 band and select a rectangular area of interest as shown in the first figure below. Then I copy the vertexes coordinates (using the “WKT-from-geometry” tool) into the “wkt” string variable of the Python code below (Code 1). This code reads the same product opened in SNAP and performs a sub-selection starting from the Polygon defined by the WKT. When looking at the size of the image obtained from the sub-selection via Snappy, I obtain a rectangle of size 20x23 (pixels) instead of 19x22 (pixels) as I have chosen in SNAP… In short, there is an unwanted addition of a column on the right and a row at the bottom (see the blu area Vs the grey-scale one of the pictures below).
This “artificial addition” seem to appear even if I build a Polygon directly in snappy (Code2). For example, if I select a 10x10 square area at the bottom-right corner of the raster (B4 rasters are 10980x10980 pixels size) when I plot the selected image it appears with an additional column.

I wonder if there is some trivial mistake in my code (I’m not expert in snappy and I have made this code mixing codes from other users and interpreting the API documentation) or if I’m doing a wrong use of the CRS of the product, or if this is simply a well-known artifact that can be solved by rounding the values at some point. I post the output of the Lat-Lon and X-Y coordinates obtained from each case study.
Below I have attached the Python codes and their outputs and some custom functions. It can be seen that the X-Y coordinates are rounded during the conversion from the lat-lot coordinates.
Other usefull informations:

  1. the area I’m studiyng is extracted from the product: “S2A_MSIL2A_20210729T100031_N0301_R122_T33TUL_20210801T154341.SAFE”.
  2. Snappy build is in it last version. I’m working on Windows 10.
  3. The geocoding of the image is the default WGS-84.


product = ProductIO.readProduct(xml_path)
B4_Product = Product('Band_04', 'Single_band', 10980, 10980)
ProductUtils.copyBand('B4', product, B4_Product, True)      
ProductUtils.copyGeoCoding(product, B4_Product)

wkt='POLYGON ((13.555187542063273 45.77662734455046, 13.557604620550356 45.77665788956839,'+\
' 13.557655389322349 45.774689044751646, 13.555238395838408 45.77465850181955, 13.555187542063273 45.77662734455046))'

Sub_Prod=  Subsetting(wkt, B4_Product) 
Sub_B4 = Sub_Prod.getBand('B4')    
Width = Sub_B4.getRasterWidth()
Height = Sub_B4.getRasterHeight()
print('Sub Selection size = ',Width,Height)

Sub_B4_data = np.zeros( (Height, Width), dtype = np.float32)
plt.figure(figsize = (6,6))
plt.imshow(Sub_B4_data, cmap = cm.gray)

B4gc = B4_Product.getSceneGeoCoding()
LatLon = wkt[10:-3].split(' ')

for lines in range(0,len(LatLon)-2,2):
        GeoLatLon = GeoPos(float(LatLon[lines+1][:-1]),float(LatLon[lines]))
        Pixel_XY = B4gc.getPixelPos(GeoLatLon, None)
        print(str(lines)+':', 'Pixel(', Pixel_XY.getX(),',', Pixel_XY.getY(),')', \
              '  LatLon(',GeoLatLon.getLon(), ',',GeoLatLon.getLat(), ')')


XY_Pairs = [PixelPos(10970,10970),PixelPos(10980,10970),PixelPos(10980,10980),PixelPos(10970,10980)]

wkt, LatLon = WKT_PixelToGeo(XY_Pairs,B4_Product)


XY_Back = []
for i in range(len(LatLon)): 
    Pixel = XY_From_LatLon(LatLon[i], B4_Product)
    print(str(i), ': ', 'Pixel(', Pixel.getX(),',', Pixel.getY(),')', LatLon[i].getLon(), LatLon[i].getLat())

sub_product= Subsetting(wkt, B4_Product)
sub_B4 = sub_product.getBand('B4')    
width = sub_B4.getRasterWidth()ttr
height = sub_B4.getRasterHeight()
print('Sub Selection size = ',width,height)
sub_B4_data = np.zeros(width*height,dtype = np.float32)
sub_B4_data.shape = height, width
fig = plt.imshow(sub_B4_data, cmap = cm.gray)

Custom functions:

def XY_From_LatLon(LatLon, prod):
    Pgc =  prod.getSceneGeoCoding()
    Pixel = Pgc.getPixelPos(LatLon, None)
    return Pixel

def WKT_PixelToGeo( XY_Polygon_Vces, prod):
    LatLonStr = []
    LatLon = []
    Pgc = prod.getSceneGeoCoding()
    for i in range(len(XY_Polygon_Vces)):
        Geopos=  Pgc.getGeoPos(XY_Polygon_Vces[i], None)  
        LatLonStr.append(str(Geopos.getLon())+' '+str(Geopos.getLat()))
        print(str(i), ': ', 'Pixel(', XY_Polygon_Vces[i].getX(),',',XY_Polygon_Vces[i].getY(),')', Geopos.getLon(), Geopos.getLat())
    chain = ', '.join(LatLonStr)+ ', '+ LatLonStr[0]
    WKT = 'POLYGON (('+ chain + '))'

    return WKT, LatLon

def Subsetting(wkt, prod):
    geom = WKTReader().read(wkt)
    op = SubsetOp()
    Sub_Product = op.getTargetProduct() 
    return Sub_Product

Outputs from Code1:

Outputs from Code2:

Thanks for the exhaustive documentation of the issue. That’s very appreciated.

The issue is not in your implementation but in SNAP and should be fixed with the next major release.
The issue description is about multi-size but when looking at the solution also single size are affected.

You could add another subset operation to overcome the problem in the short-term.
There you could define the pixel region. Retrieve the width and height from the previous subset result and subtract one.
This should give you the desired result.

Thanks @marpet for the response, I have made the solution you suggested and everything is ok.