Iterating over shapefile to mask and extract raster data in Python

Question

I have a shapefile that is a grid with multiple polygons that have a specific id in them ranging from 1 to 50. Each polygon overlays an agriculture research plot that I would like to calculate vegetation indices. I then want to append this information to an Excel file.

But I am having trouble with getting the loop to work, where only the first polygon appends data to the Excel sheet, it seems like it doesn't even loop through. Any ideas or improvements?

import rasterio
import fiona
import numpy as np
import shapefile as shp
from openpyxl import load_workbook

with fiona.open("path/to/shapefile.shp", "r") as shapefile:
    shpfl = shp.Reader("path/to/shapefile.shp")
    recs = shpfl.shapeRecords()
    count = 0
    wb = load_workbook(wb_output)
    ws = wb.active
    for feature in shapefile:
        shape =[feature["geometry"]]
        plotNum = recs[count].record[2]
        sheetRef = plotNum + 1

        #open and mask the multispectral image to the plot
        with rasterio.open(multispec) as src:
            out_image, out_transform = rasterio.mask.mask(src, shape, crop = True)
            out_meta = src.meta
            
            #extract color bands
            red = out_image[0,:,:].astype('float64')
            nir = out_image[3,:,:].astype('float64')
            red_edge = out_image[4,:,:].astype("float64")
                
            #ndvi           
            num = np.subtract(nir, red)
            denom = np.add(nir, red)
            np.seterr('ignore', 'ignore')
            ndvi = np.divide(num, denom)
            
            flatndvi = ndvi.flatten()
            avg_ndvi = np.nanmean(flatndvi)
            

            cell = ws['c' + str(sheetRef)]
            cell.value = avg_ndvi
            wb.save(wb_output)

Answer 1

As @bugmenot123 stated, no need to use shapefile module and Fiona has all you need.

import rasterio
import fiona
import numpy as np
from openpyxl import load_workbook
from rasterio.mask import mask

wb_output = "path/to/excelfiel.xlsx"
multispec = "path/to/image.tif"
shp_file = "path/to/shapefile.shp"
plot_num_field_name = "PLOTNUM_FIELDNAME_in_SHAPEFILE"

with fiona.open(shp_file, "r") as shapefile:
    wb = load_workbook(wb_output)
    ws = wb.active

    for feature in shapefile:
        shape = [feature["geometry"]]
        sheetRef = feature["properties"][plot_num_field_name] + 1

        #open and mask the multispectral image to the plot
        with rasterio.open(multispec) as ras:
            out_image, out_transform = mask(ras, shape, crop=True)

            #extract color bands
            red = out_image[0,:,:].astype('float64')
            nir = out_image[3,:,:].astype('float64')
            red_edge = out_image[4,:,:].astype("float64")

            #ndvi
            num = np.subtract(nir, red)
            denom = np.add(nir, red)
            np.seterr('ignore', 'ignore')
            ndvi = np.divide(num, denom)

            flatndvi = ndvi.flatten()
            avg_ndvi = np.nanmean(flatndvi)

            cell = ws['c' + str(sheetRef)]
            cell.value = avg_ndvi
            wb.save(wb_output)

Toy data result: