title: How to calculate average pixel values for multiple bands that fall within polygons
body: <p>I have a hyperspectral raster image of an agricultural field with 270 spectral bands. I created a polygon shapefile that delineates which pixels belong to each treatment. There are 250 individual polygons that each correspond to a replication of each treatment. I want to find the average pixel value for each band for all of the pixels that fall within each polygon.</p>
  
  <p><a href="https://i.stack.imgur.com/tbTcP.jpg" rel="nofollow noreferrer">Image of raw hyperspectral data</a></p>
  
  <p><a href="https://i.stack.imgur.com/jo8EE.jpg" rel="nofollow noreferrer">Image of polygons delineating treatments</a></p>
  
  <p>I tried using the zonal statistics tool in both ArcGIS and QGIS but both tools can only run statistics on one band at a time. Doing this 270 times seems a little excessive.</p>
  
  <p>I also tried to use the split raster tool in ArcGIS to divide the raster into 250 individual rasters corresponding to each polygon. Once I split the raster, I tried using the band collection statistics tool but found that I could not input all rasters simultaneously although the tool appears to be capable of doing so. Each attempt results in the following error: ERROR 000964  Specified extent is invalid.</p>
  
  <p>I've been conducting my analyses in ArcGIS Pro, QGIS (v.3.4.11), and Python (v.3.7) primarily using GDAL. So, I am open to using any of these options to conduct further analysis. I think this might be doable in Python, but my coding skills aren't great and I'm not sure where to start.</p>