How to extract pixel value of every band in a rasterstack using a polygon?

Question

I am trying to extract the pixel values of a raster from a polygon that has several features. My raster is a rasterstack with 4 bands. I have found how to do it for a whole raster, but I need the info PER BAND. Any hints?

from rasterstats import zonal_stats
import os
import numpy as np 
import statistics
 
shapefile = Class1.shp
geotiff = tile_105.tif

# calculate all zonal stats
stats = zonal_stats(shapefile, geotiff)

# extract the mean and store it
single_mean = [f['mean'] for f in stats]

val_list = []

# only store the positive values. 
for val in single_mean:
    if val != None :
        val = [float(val)]
      
        val_list.append(val[0])

all_mean = statistics.mean(val_list)

Answer 1

rasterio to read a [multi-band raster]!

https://rasterio.readthedocs.io/en/latest/topics/reading.html

import rasterio

dataset = rasterio.open('path_to_your_multiband.tif')

You may also want to define the extent of your AOI by defining the cropping window ((row_start, row_stop), (col_start, col_stop))

window = ((20, 50), (10, 40))
val_list[]
with rasterio.open(dataset) as src:
    # Create zero array (you may want to set dtype too) for negative values
    # only store the positive values. 
    array = np.zeros((window[0][1] - window[0][0],
                      window[1][1] - window[1][0],
                      len(bands)))
    # Fill the array
    for i, band in enumerate(bands):
        array[:,:,i] = src.read(band, window=window) 
        val_list.append(i[0])

all_mean = statistics.mean(val_list)