topic Re: Loop in python with Zonal Statistics in Python Questions

Loop in python with Zonal Statistics

KrzysztofSterenczak — Thu, 13 Mar 2014 12:58:09 GMT

Hi,

I have following problem. I have a layer with thousands of objects - crowns. I have CHM (nDSM) with "height" values. I would like to calculate "All" height statistics for crowns based on CHM, crown after a crown - because to do so directly in Spatial Analyst is impossible (to many objects in layer). I tried to write a Python script (one of my first ones) but it simply does not work:

import arcpy
import os
from arcpy import env
from arcpy.sa import *
# Set environment settings
env.workspace = "F:\\_IBL\\Projekty\\GDLP_Sudety_Beskidy\\Analizy\\ALS2007vs2012\\analiza_probna\\baza.gdb"
outWorkspace = "F:\\_IBL\\Projekty\\GDLP_Sudety_Beskidy\\Analizy\\ALS2007vs2012\\analiza_probna\\baza.gdb"
featureClass = "korony_proba" # crowns
raster2007 = "T2007_proba" #CHM (nDSM) ??? 32 bit, flout

rows = arcpy.SearchCursor(featureClass)
row = rows.next()
for row in rows:
inZoneData = "featureClass"
zoneField = "ID_OK"
inValueRaster = "raster2007"
outTable = "stat_2007"
outZSaT = ZonalStatisticsAsTable(inZoneData, zoneField, inValueRaster, outTable, "NODATA", "ALL")

I have following error on the beginning:
Runtime error Traceback (most recent call last):
File "<string>", line 11, in <module>
File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\sa\Functions.py", line 6195, in ZonalStatisticsAsTable     statistics_type)
File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\sa\Utils.py", line 47, in swapper     result = wrapper(*args, **kwargs)
File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\sa\Functions.py", line 6187, in Wrapper     statistics_type)
File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\geoprocessing\_base.py", line 498, in <lambda>     return lambda *args: val(*gp_fixargs(args, True)) ExecuteError: ERROR 000865: Input raster or feature zone data: featureClass does not exist. ERROR 001000: Zone field: Field ID_OK does not exist ERROR 000865: Input value raster: raster2007 does not exist.
Can you help me please.
Thanks,
Kris

Re: Loop in python with Zonal Statistics

ClintDow — Thu, 13 Mar 2014 21:22:06 GMT

Hello, a couple things that might help:
You are mixing up two methods of using the cursor in this section

rows = arcpy.SearchCursor(featureClass) row = rows.next() for row in rows:

You do not need to assign row = rows.next() unless you are using a while loop, in a for loop using the cursor as a sequence to loop through is enough.

In the line inZoneData = "featureClass" you are assigning the string 'featureClass' rather than assigning the value from the featureClass variable. Therefore the script is looking for a feature class named 'featureClass' in your workspace rather than a feature class named 'korony_proba' - this is the cause of ERROR 000865: Input raster or feature zone data: featureClass does not exist.

Since that doesn't exist, there is no ID_OK field to find, hence the ERROR 001000: Zone field: Field ID_OK does not exist

ERROR 000865: Input value raster: same issue as with the feature class, its looking for a raster named 'raster2007' rather than assigning 'T2007_proba' to that variable.

Re: Loop in python with Zonal Statistics

KrzysztofSterenczak — Fri, 14 Mar 2014 07:37:38 GMT

Dear Clinton,
Thank you very much for tis tips. Now I solved a problem.
Final code:

import arcpy
import os
from arcpy import env
from arcpy.sa import *

# Set environment settings
env.workspace = "F:\\_IBL\\Projekty\\GDLP_Sudety_Beskidy\\Analizy\\ALS2007vs2012\\analiza_probna\\baza.gdb"
outWorkspace = "F:\\_IBL\\Projekty\\GDLP_Sudety_Beskidy\\Analizy\\ALS2007vs2012\\analiza_probna\\baza.gdb"

rows = arcpy.SearchCursor(featureClass)
for row in rows:
# Set local variables
inZoneData = "rosnace"
zoneField = "ID_OK"
inValueRaster = "T2007"
outTable = "stat_2007_test1"
# Execute ZonalStatisticsAsTable
outZSaT = ZonalStatisticsAsTable(inZoneData, zoneField, inValueRaster, outTable, "NODATA", "ALL")

As a result I have table with all statistics and with error:
ERROR 000872: Output table: Output F:\_IBL\Projekty\GDLP_Sudety_Beskidy\Analizy\ALS2007vs2012\analiza_probna\baza.gdb\stat_2007_test1 exists. It cannot be overwritten since overwrite is off.

When I added:
arcpy.env.overwriteOutput = True
Error did not show off but table was empty and procedure did not stoped.
Nevertheless I reached my aim so thank you very much ones again.
Kris

Re: Loop in python with Zonal Statistics

XanderBakker — Fri, 10 Dec 2021 22:36:02 GMT

Hi Krzysztof,

Clint already indicated some flaws in the code. Another one is the fact that you do a loop over your features, but nothing is actually done with the feature itself. What you should do in the loop is create a feature layer of the single feature and use that feature layer in the zonal statistics tool.

See below an example of how this could be done:

import arcpy
import os

# Set environment settings
outWorkspace = r"C:\Project\_Forums\_json\fgdb\crawlAHN2.gdb"
arcpy.env.workspace = outWorkspace
arcpy.env.overwriteOutput = True

# reference to the data
fc = os.path.join(outWorkspace, "crowns")
raster = os.path.join(outWorkspace, "AHN2r")
zoneField = "ID_OK"

# temporal output ZS table will be (over)written to in memory workspace
outTable = os.path.join("IN_MEMORY", "stat_zs")

# name of feature layer (lives in memory)
flyr_name = 'a_crown'

# get list of zones (this is used for the loop)
lst_zones = [row[0] for row in arcpy.da.SearchCursor(fc, (zoneField))]

# create an empty dictionary to store the statistics results
dct_stats = {}

# create a dictionary that holds the statistics type and the corresponding output field names
dct_stat_fld = {'MEAN': 'cr_mean',
                'MAX':  'cr_max',
                'MIN': 'cr_min',
                'STD': 'cr_std'}

# list of fields in ZS table that hold the statistics I'm interested in ['MEAN','MIN','MAX','STD']
lst_flds = [fld for fld in dct_stat_fld.iterkeys()]

# get sa license
arcpy.CheckOutExtension("Spatial")

# now loop through zones and make a featurelayer with only that zone feature (crown)
cnt = 0
for zone in lst_zones:
    cnt += 1
    # give some feedback on progress
    if cnt % 10 == 0:
        print "Processing crown: {0} ({1})".format(cnt, zone)

    # create the where clause for that crown
    where = "{0} = '{1}'".format(arcpy.AddFieldDelimiters(fc, zoneField), zone)

    # make the featurelayer of current crown (zone)
    arcpy.MakeFeatureLayer_management(fc, flyr_name, where)

    # calculate the zonal statistics for that crown
    arcpy.sa.ZonalStatisticsAsTable(flyr_name, zoneField, raster, outTable, "NODATA", "ALL")

    # read the ZS table and write to dictionary
    with arcpy.da.SearchCursor(outTable, lst_flds) as curs:
        for row in curs:
            i = 0
            for fld in lst_flds:
                key = (zone, fld) # looks like ('crown ID', 'statistics type like MIN')
                dct_stats[key] = row
                i += 1

# return sa license
arcpy.CheckInExtension("Spatial")

# create update field collection (field  for the update cursor)
flds = [zoneField]
flds.extend([dct_stat_fld[fld] for fld in lst_flds])

# do an update cursor to write the results from dictionary to crown featureclass
with arcpy.da.UpdateCursor(fc, flds) as curs:
    for row in curs:
        zone = row[0]
        i = 1
        for fld in lst_flds:
            key = (zone, fld)
            if key in dct_stats:
                row = dct_stats[key]
            i += 1
        curs.updateRow(row)

I hope this is useful for you.

Kind regards,

Xander