How to calculate average pixel values for multiple bands that fall within polygons

Question

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.

Image of raw hyperspectral data

Image of polygons delineating treatments

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.

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.

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.

I don't have a direct answer for you, but perhaps this link will help you to get more understanding of coding with QGIS. Good luck! — Johan

fruitspunchsamurai fruitspunchsamurai · Accepted Answer · 2019-10-11T09:37:40

You can use Spectral Python package to read your hyperspectral image. Lets say,

from spectral import*
img=open_image('<file location>')  #read into img

make a numpy mask image from the polygon shapefile you created, probably according to the classes(color coding).From this image, you can extract the pixel co-ordinates.

import numpy as np
from matplotlib.image import imread
mask_img=imread('<polygon mask image location>')

lets say for one polygon you have assigned class 1,so for all the pixel points in class1 polygon

x,y=np.where(mask_img==1)  #get the required co-ordinates

use the same on the hyperspectral image

band_avg=[]
for band in np.arange(270):
    sum,avg=0.0,0.0
    for i in np.arange(len(x)):   #no. of pixels in polygon class1
        sum+= img[x[i],y[i]]
     avg=sum/len(x)
     band_avg.append(avg)
print(band_avg)

band_avg will return average pixel value for each band for all of the pixels in a particular class/polygon color. You can repeat the above for other classes/polygons by getting x,y for different class ID's.

How to calculate average pixel values for multiple bands that fall within polygons

1 Answers