After import a raster file and add buffer of 20 meters I want to write a new raster with this dimensions, but getting error:
import rasterio
from rasterio import features
from rasterio import transform
from shapely.geometry import shape, Point, Polygon
import numpy as np
import zipfile,contextlib
from io import BytesIO
from contextlib import ExitStack
#read cloud raster from zip
with zipfile.ZipFile('...img/S2B_MSIL2A_20230926T132239_N0509_R038_T22KHC_20230926T171115.SAFE.zip',"r") as zf_:
zfiledata=BytesIO(zf_.read('S2B_MSIL2A_20230926T132239_N0509_R038_T22KHC_20230926T171115.SAFE/GRANULE/L2A_T22KHC_A034241_20230926T132553/IMG_DATA/R20m/T22KHC_20230926T132239_SCL_20m.jp2')) with contextlib.ExitStack() as stack:
uzFile=stack.enter_context(zf.open('S2B_MSIL2A_20230926T132239_N0509_R038_T22KHC_20230926T171115.SAFE/GRANULE/L2A_T22KHC_A034241_20230926T132553/IMG_DATA/R20m/T22KHC_20230926T132239_SCL_20m.jp2'))
_uzFileBytes=stack.enter_context(BytesIO(_uzFile.read()))
with rasterio.open(_uzFileBytes) as src_raster:
raster_pixels=src_raster.read(1)
raster_transform=src_raster.transform
raster_meta=src_raster.meta scl_pixels=raster_pixels.copy()
#detect pixels with cloud and shadow
index_pixels=np.where(np.isin(scl_pixels.astype('float32'),[3,8,9,10]),scl_pixels.astype('float32'),None)
cldsdw_rows,cldsdw_cols=np.where(np.array(index_pixels,dtype=np.float32)>0)
cldsdw_coords = raster_transform(cldsdw_rows,cldsdw_cols)
cldsdw_points = [Point(x, y) for x, y in zip(cldsdw_coords[0], cldsdw_coords[1])]
cldsdw_polygon = Polygon(cldsdw_points)
cldsdw_bufygon = cldsdw_polygon.buffer(20)
##add buffer 20m
cldsdw_features = [{'geometry': shape(geom)} for geom in cldsdw_bufygon.geoms]
with rasterio.open('...img/r20moutput.tif', 'w', *raster_meta) as dst:
rasterio.features.rasterize(cldsdw_features, dst.write, transform=raster_transform)
but on write, its getting errors like TypeError:
... /opt/.virtualenv/visaopq/lib/python3.12/site-packages/rasterio/features.py:336: ShapeSkipWarning: Invalid or empty shape {'geometry': <POLYGON ((799973.998 7735779.078, 800020 7735793.552, 800020 7735840, 80002...>} at index 0 will not be rasterized. warnings.warn( /opt/.virtualenv/visaopq/lib/python3.12/site-packages/rasterio/features.py:336: ShapeSkipWarning: Invalid or empty shape ...warnings.warn( /opt/.virtualenv/visaopq/lib/python3.12/site-packages/rasterio/features.py:336: ShapeSkipWarning: Invalid or empty shape {'geometry': <POLYGON ((903779.055 7780200.043, 903777.76 7780200.126, 903775.841 7780200...>} at index 22 will not be rasterized. warnings.warn( Traceback (most recent call last): File "", line 2, in File "/opt/.virtualenv/visaopq/lib/python3.12/site-packages/rasterio/env.py", line 410, in wrapper return f(*args, **kwds) ^^^^^^^^^^^^^^^^ File "/opt/.virtualenv/visaopq/lib/python3.12/site-packages/rasterio/features.py", line 378, in rasterize if len(out_shape) != 2: ^^^^^^^^^^^^^^ TypeError: object of type 'method' has no len()
some outputs to clarify:
scl_pixels
array([[ 4, 4, 4, ..., 5, 5, 5],
[ 4, 4, 4, ..., 5, 5, 5],
[ 4, 4, 4, ..., 5, 5, 5],
...,
[10, 10, 10, ..., 4, 4, 4],
[10, 10, 10, ..., 4, 4, 4],
[10, 10, 10, ..., 4, 4, 4]], dtype=uint8)
index_pixels
array([[None, None, None, ..., None, None, None],
[None, None, None, ..., None, None, None],
[None, None, None, ..., None, None, None],
...,
[10.0, 10.0, 10.0, ..., None, None, None],
[10.0, 10.0, 10.0, ..., None, None, None],
[10.0, 10.0, 10.0, ..., None, None, None]], dtype=object)
cldsdw_features
[{'geometry': <POLYGON ((799973.998 7735779.078, 800020 7735793.552, 800020 7735840, 80002...>}, {'geometry': <POLYGON ((870280.908 7799934.063, 870280.861 7799934.194, 870280.384 779993...>}, {'geometry': <POLYGON ((870329.993 7799782.7, 870328.655 7799783.529, 870327.117 7799784....>}, {'geometry': <POLYGON ((870470.013 7799982.7, 870469.057 7799983.259, 870467.453 7799984....>}, {'geometry': <POLYGON ((871560 7758302.129, 871602.469 7758315.491, 871600 7752439.992, 8...>}, {'geometry': <POLYGON ((874818.993 7799059.131, 874780.026 7799818.976, 874780 7799820, 8...>}, {'geometry': <POLYGON ((875072.266 7799664.207, 875070.737 7799663.127, 875069.05 7799662...>}, {'geometry': <POLYGON ((875580 7800015.404, 875580 7800040, 875580.247 7800040, 875600 78...>}, {'geometry': <POLYGON ((879300 7760737.359, 879351.351 7760753.516, 879340 7711819.995, 8...>}, {'geometry': <POLYGON ((886740 7763078.201, 886790.736 7763094.164, 886780 7720279.995, 8...>}, {'geometry': <POLYGON ((887160 7763210.345, 887211.078 7763226.416, 887200 7717519.995, 8...>}, {'geometry': <POLYGON ((901103.816 7767597.478, 901200 7767627.74, 901200 7790360, 901200...>}, {'geometry': <POLYGON ((901299.948 7790439.169, 901282.326 7790580.148, 901282.042 779058...>}, {'geometry': <POLYGON ((901320 7790380.018, 901320 7790480, 901320.247 7790480, 901340 77...>}, {'geometry': <POLYGON ((901700.483 7789500.028, 901699.744 7789500.002, 901697.801 778950...>}, {'geometry': <POLYGON ((901845.698 7789420.85, 901840.003 7789779.683, 901840.249 7789779...>}, {'geometry': <POLYGON ((901890.001 7788502.704, 901888.889 7788503.371, 901887.312 778850...>}, {'geometry': <POLYGON ((901959.08 7789720.045, 901958.04 7789720.096, 901956.098 7789720....>}, {'geometry': <POLYGON ((902082.445 7775980.168, 902081.96 7775980.096, 902080 7775980, 90...>}, {'geometry': <POLYGON ((902259.289 7788783.873, 902224.383 7788827.506, 902223.253 778882...>}, {'geometry': <POLYGON ((902477.653 7791640.143, 902477.556 7791640.15, 902475.654 7791640...>}, {'geometry': <POLYGON ((902849.006 7768146.566, 903497.172 7768350.498, 903239.994 775753...>}, {'geometry': <POLYGON ((903779.055 7780200.043, 903777.76 7780200.126, 903775.841 7780200...>}]
cldsdw_bufygon
<MULTIPOLYGON (((799973.998 7735779.078, 800020 7735793.552, 800020 7735840,...>
Im trying to convert cldsdw_bufygon to pixels like on index_pixels and write an geotif image. Or overlay the multipolygon on original image...
