I have a shapefile of thousands of polygons. I'm looking to loop through each polygon and create a raster of the shapefiles boundaries for each polygon. I'm looking for each raster to have the value of 1 for the areas where the polygon exists and a value of 0 where it doesn't. I've tried the below but the rasters are off both in terms of shape and crs.
eco_sys = gpd.read_file(all_habitats_path)
# eco_sys = eco_sys[eco_sys.EUNIScombD == 'A6.51']
eco_len = len(eco_sys.geometry)
eco_list = [1] * eco_len
rst = rasterio.open(cropped_rast_template)
transform, width, height = calculate_default_transform(
rst.crs, ENV.get("jose_crs_str"), rst.width, rst.height, *rst.bounds
)
meta = rst.meta.copy()
meta.update({
'compress': 'lzw',
'crs': ENV.get("jose_crs_str"),
'transform': transform,
'width': width,
'height': height,
'nodata': -999
})
# eco_sys.geometry.values[0].buffer(0)
shapes = ((geom, value)
for geom, value in zip(eco_sys.geometry, eco_list))
with rasterio.open('data/tmp/all_habs.tif', 'w+', **meta) as out:
out_arr = out.read(1)
# combined_rasters = None
for idx, shape in enumerate(shapes):
print(' idx, shape : ', idx, shape)
if idx != 11:
burned = features.rasterize(
shapes=shape, fill=0, default_value=1, out=out_arr, transform=out.transform)
with rasterio.open('data/tmp/eco_sys/'+str(idx)+'.tif', 'w+', **meta) as dst:
reproject(
source=burned,
destination=rasterio.band(dst, 1),
src_transform=rst.transform,
src_crs=rst.crs,
dst_transform=transform,
dst_crs=ENV.get("jose_crs_str"),
resampling=Resampling.nearest,
dst_nodata=-999)
dst.write(burned, 1)