Hi,
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:
- the area I’m studiyng is extracted from the product: “S2A_MSIL2A_20210729T100031_N0301_R122_T33TUL_20210801T154341.SAFE”.
- Snappy build is in it last version. I’m working on Windows 10.
- The geocoding of the image is the default WGS-84.
Code1
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)
Sub_B4.readPixels(0,0,Width,Height,Sub_B4_data)
plt.figure(figsize = (6,6))
plt.imshow(Sub_B4_data, cmap = cm.gray)
plt.show()
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(), ')')
Code2
XY_Pairs = [PixelPos(10970,10970),PixelPos(10980,10970),PixelPos(10980,10980),PixelPos(10970,10980)]
wkt, LatLon = WKT_PixelToGeo(XY_Pairs,B4_Product)
print('-'*40)
XY_Back = []
for i in range(len(LatLon)):
Pixel = XY_From_LatLon(LatLon[i], B4_Product)
XY_Back.append(Pixel)
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.readPixels(0,0,width,height,sub_B4_data)
sub_B4_data.shape = height, width
plt.figure('Test_Rectangle')
fig = plt.imshow(sub_B4_data, cmap = cm.gray)
fig.axes.get_xaxis().set_visible(True)
fig.axes.get_yaxis().set_visible(True)
plt.show()
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)
LatLon.append(Geopos)
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()
op.setSourceProduct(prod)
op.setGeoRegion(geom)
Sub_Product = op.getTargetProduct()
return Sub_Product
Outputs from Code1:
Outputs from Code2:
