Python Zonal Statistics – Generating Zonal Stats for Multiple Rasters Using Rasterstats

pythonrasterstatszonal statistics

I'm trying to generate the zonal stats for multiple rasters without using ArcPy.

I've successfully obtained results using rasterstats for one image and one polygon. I would like to use the same polygon but many rasters. My code is below and not working. Does anyone know what I'm doing wrong?

Lista_raster = [
   "NDVI_2019-05-05.tiff", "NDVI_2019-05-25.tiff", "NDVI_2019-06-04.tiff", 
   "NDVI_2019-06-09.tiff", "NDVI_2019-06-14.tiff", "NDVI_2019-06-19.tiff",
   "NDVI_2019-06-24.tiff", "NDV_2019-05-05.tiff"
]
for raster in Lista_raster:
    stats = zonal_stats("Pivo_4_wgs.shp", Lista_raster,
            stats=['min', 'max', 'median', 'majority', 'std'],  
            geojson_out=True)
result = {"type": "FeatureCollection","features": stats}
import json
with open('lista.geojson', 'w') as outfile:
   json.dump(result, outfile) 

Best Answer

Notice this part of your code:

for raster in Lista_raster:
    stats = zonal_stats("Pivo_4_wgs.shp", Lista_raster,  # <---here

You're passing the entire list for each iteration of your loop. Try passing the the element of the iterator instead:

for raster in Lista_raster:
    stats = zonal_stats("Pivo_4_wgs.shp", raster, ...

You're also only writing a single file at the end of the loop, instead of a file for each result. So your code becomes:

import json
from pathlib import Path

Lista_raster = [
   "NDVI_2019-05-05.tiff", "NDVI_2019-05-25.tiff", "NDVI_2019-06-04.tiff", 
   "NDVI_2019-06-09.tiff", "NDVI_2019-06-14.tiff", "NDVI_2019-06-19.tiff",
   "NDVI_2019-06-24.tiff", "NDV_2019-05-05.tiff"
]

for raster in Lista_raster:
    stats = zonal_stats("Pivo_4_wgs.shp", raster,
            stats=['min', 'max', 'median', 'majority', 'std'],  
            geojson_out=True)
    result = {"type": "FeatureCollection", "features": stats}

    outname = Path(raster).name + '.geojson'
    with open(outname, 'w') as outfile:
        json.dump(result, outfile)