I have a Gtiff layer and a shapefile with 5 polygons that are placed across the Gtiff (both have the same crs). I'm trying to create a NumPy array for each polygon that contains the pixels that are masked by that polygon. The code I managed to build is attaching the entire raster to each polygon and not just the relevant pixels. Any suggestions? polygons = gpd.read_file(poly_path.shp) burn_attribute = "id" # should be unique for the polygons and not contain zero with rasterio.open(raster_path.tif") as src: raster = src.read(1) pixel_values = [] for i in range(len(polygons)): # Create a mask from the current polygon mask = rasterio.features.geometry_mask([polygons.geometry], out_shape=raster.shape, transform=src.transform, all_touched=True, invert=False) # Extract the pixel values within the polygon polygon_values = raster[mask] # Append the pixel values to the list pixel_values.append(polygon_values) # Add the pixel values as a new column to the polygons GeoDataFrame polygons['pixel_values'] = pixel_values extracted_values = raster[polygons] Continue reading...