Python Zonal Statistics Problems in QGIS

Question

I am trying to run Zonal Stats for Block Group shape files against NLCD Percent Impervious raster. I have shapefiles for 48 contiguous states all projected to the same as the raster.

I've written this Python script to try and do each state shapefile in a loop. At a Linux command line I run the script, it prints out all the shapefile names to the console and does not an error but the Zonal Stats are never added to the shapefiles. It lists all the shapefile names in a second.

If I try running it from the QGIS Python Console it appears to run the script using execfile() and no error or print statements are given but the Zonal Stats are also not added to the shapefiles. Does anyone know what I am doing wrong? Zonal Stats is running fine from the QGIS desktop.

I see there are a number of other topics on this subject but I haven't found a response in them that solves this issue.

import os, sys, traceback
from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.core import *
from qgis.utils import iface
from qgis.analysis import QgsZonalStatistics

def main():

    try:

        shp_fldr = "/media/psf/Home/Documents/Projects/NLCD/blocks_albers"

        raster = "/media/psf/Home/Documents/Projects/data/2016/NLCD_2016_Impervious_L48_20190405.img"
        f_lst = os.listdir(shp_fldr)

        for fl in f_lst:

            if fl.find(".shp") > 0 and fl.find(".xml") < 0:

                print "start: " + fl

                shp_fl = os.path.join(shp_fldr,fl)

                #specify polygon shapefile vector

                polygonLayer = QgsVectorLayer(shp_fl, fl, "ogr")

                polygonLayer.startEditing()

                # usage - QgsZonalStatistics (QgsVectorLayer *polygonLayer, const QString &rasterFile, const QString &attributePrefix="", int rasterBand=1)

                zoneStat = QgsZonalStatistics (polygonLayer, raster, 'ipv-', 1, QgsZonalStatistics.Mean)

                zoneStat.calculateStatistics(None)

                polygonLayer.commitChanges()

                print "end: " + fl

    except:

        # Get the traceback object

        tb = sys.exc_info()[2]

        tbinfo = traceback.format_tb(tb)[0]

        # Concatenate information together concerning the error into a message string

        pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1])

        # Return python error messages for use in script tool or Python Window

        print pymsg

if __name__ == '__main__':

    main()

Answer 1

I got it to work. I stripped the script down to just the basics and ran it in the QGIS Python Console and it ran through all the state Block Group shapefiles just as I'd expect. The Zonal Stats fields are now in the attribute tables.

import os
from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.core import *
from qgis.utils import iface
from qgis.analysis import QgsZonalStatistics

shp_fldr = "/media/psf/Home/Documents/Projects/NLCD/blocks_albers"
raster = "/media/psf/Home/Documents/Projects/data/NLCD/2016/NLCD_2016_Impervious_L48_20190405.img"

f_lst = os.listdir(shp_fldr)

for fl in f_lst:

    if fl.find(".shp") > 0 and fl.find(".xml") < 0:

        print "start: " + fl
        shp_fl = os.path.join(shp_fldr,fl)

        polygonLayer = QgsVectorLayer(shp_fl, fl, "ogr")
        polygonLayer.startEditing()
        zoneStat = QgsZonalStatistics (polygonLayer, raster, 'ipv-', 1, QgsZonalStatistics.Mean)

        zoneStat.calculateStatistics(None)
        polygonLayer.commitChanges()

        print "end: " + fl

Answer 2

In your code you just create layer. To visualize it in qgis map, you have to add the layer to the map canvas. https://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/loadlayer.html

QgsProject.instance().addMapLayer(layer)

Or

iface.addVectorLayer...