3

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)
2
  • 1
    You are mixing the fiona and shapefile packages for no reason. Fiona has all you need, so try remove shapefile entirely. Then it will be easier to see what is happening. – bugmenot123 Jan 16 at 8:04
  • For debugging in general it would be good to either focus on the iteration through the polygons and their data extraction OR the writing of data into an excel sheet. Right now both things might be broken. If you print avg_ndvi each time, does it update? Does sheetRef? – bugmenot123 Jan 16 at 8:08
4

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:

enter image description here

1
  • thanks for the help. it worked well! – Mlien Jan 18 at 0:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.