Batch Zonal Statistics Script - Error 010423

[N__R_Wilson]

I created a script to batch process Zonal Statistics as Table with user inputs, option to add a field, and a single appended table as an output.

It runs fine with some data. But with other data it throws a ERROR 010423: N:\Gabions\NDVI\35_38_May_June\t_t324 does not have valid statistics as required by the operation
When I modify the code to run from PyScripter it works fine and doesn't throw the error. (Modifications: hard-code the inputs, change AddMessage/AddWarning to print.) I haven't been able to find much information on this error within a python setting.

I tried adding Calculate Statistics just before the Zonal Statistics as Table but nothing changed...
arcpy.CalculateStatistics_management(raster,"","","","OVERWRITE",zonefile)

I'm pretty new to scripting. Any help?

# Import modules
import arcpy
import os
import shutil

# Get parameters
indir = arcpy.GetParameterAsText(0)
zonefile = arcpy.GetParameterAsText(1)
addfield = arcpy.GetParameterAsText(2)
startcharin = arcpy.GetParameterAsText(3)
endcharin = arcpy.GetParameterAsText(4)
outdir = arcpy.GetParameterAsText(5)
outfile = arcpy.GetParameterAsText(6)

# Set variables for zonal statistics and add field routines
zonefield = "UID"
startchar = int(startcharin) - 1
endchar = int(endcharin) - 1
outtable = ""
year = ""

# Set global environments and make scratch
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True
scratchdir = "C:\\TEMP\\Scratch\\ZoneStatsBatch"
if not os.path.exists(scratchdir):
    os.makedirs(scratchdir)

# Check the input data
fieldlist = arcpy.ListFields(zonefile,"*","ALL")
fieldnames = []
for field in fieldlist:
    fieldnames.append(field.name)
if not "UID" in fieldnames:
    arcpy.AddWarning(zonefield + " field does not exist in " + zonefile)
arcpy.env.workspace = indir
rasters = arcpy.ListRasters("*")
if not rasters:
    arcpy.AddWarning(indir +" contains no rasters \nZonal statistics failed")

# Zonal Statistics
index = 0
for raster in rasters:
    outtable = scratchdir + "\\scratchouttable" + str(index) + ".dbf"
    arcpy.sa.ZonalStatisticsAsTable(zonefile,zonefield,raster,outtable,"DATA","ALL")
    arcpy.AddMessage("Zonal statistics have been calculated for raster: " + raster)
# Add Field and Calculate Field if addfield boolean = True
    if addfield:
        try:
            year = raster[startchar:endchar]
            arcpy.AddField_management(outtable, "Year", "TEXT", "", "", 4)
            arcpy.CalculateField_management(outtable, "Year", year, "PYTHON")
            arcpy.AddMessage ("Year field (" + year + ") has been calculated for raster: " + raster)
        except:
            arcpy.AddWarning("Add year field failed")
    index = index + 1

# Set local variables and environment for Append
arcpy.env.workspace = scratchdir
outfiledbf = outfile +".dbf"
template = outtable
config_keyword = ""
outfullpath = outdir + "\\" + outfiledbf
schemaType ="NO_TEST"
fieldMappings = ""
subtype = ""

# Append
try:
    tablelist = arcpy.ListTables("*")
    arcpy.CreateTable_management(outdir, outfiledbf, template, config_keyword)
    arcpy.Append_management(tablelist, outfullpath, schemaType, fieldMappings, subtype)
    arcpy.AddMessage("Batch zonal statistics output table located at: " + outdir + "\\" + outfiledbf)
except:
    arcpy.AddWarning("Table append failed")

# Clean up your mess! (delete scratch)
try:
    shutil.rmtree(scratchdir,ignore_errors = "TRUE", onerror = "NONE")
    arcpy.AddMessage("Scratch workspace cleared")
except:
    arcpy.AddWarning("Scratch workspace not cleared: " + scratchdir)

[ChrisSnyder]

Some ideas:

1. Make sure your input zone rasters are integer format and have a raster attribute table

rasterObj = arcpy.Raster(r"C:\temp\myraster.img")
rasterObj.hasRAT

2. Don't use file names/directory names that start with a number.

[curtvprice]

I also recommend using arcpy functions to create and remove things so that arcpy is aware of those creations and deletions without having to invoke arcpy.RefreshCatalog().

[N__R_Wilson]

The zone file is a vector file which is allowable by the tool, it just means that ZonalStatisticsAsTable rasterizes it within the tool itself. And if the tool throws an exception about a raster it created in its little black box, is there anything to be done about it?

My data structure is solid with ESRI ArcGIS. All alpha starting characters, no spaces.

I did modify my code to include all arcpy creations and deletions, instead of switching to os and shutil. I'm still learning about what's lurking in the depths of arcpy so thanks for pointing me in that direction.

[ChrisSnyder]

is there anything to be done about it

Yes, instead of relying on the "auto-rasterization" of the Zonal Stats, do it yourself explicitly (run the FeatureToRaster/PolygonToRaster/WhateverToRaster tool yourself). Things can go wrong, if for example, your zone field in the vector layer has wonkey values such as Null. I always do this step myself as to catch/prevent issues like this... another reason to do it yourself is to explicitly control the cell alignment of the raster version of the zone vector layer.

My data structure is solid with ESRI ArcGIS. All alpha starting characters, no spaces

In your example above, you use a path of N:\Gabions\NDVI\35_38_May_June\t_t324. In this path you have a folder that starts with the number 3.

[N__R_Wilson]

!
Nice catch on the folder label. I obviously wasn't careful enough when examining my own data.

I did insert the FeatureToRaster today. Unfortunately I also realized that much of the data this script will be used for has overlapping polygons in the zone data, which led me down a totally separate rabbit hole where I need to separate each feature into it's own feature class and then rasterize those before I even get to the ZonalStats. It also extends the processing time by quite a bit...

Thanks for your input, I'm still wrestling with the beast.

[N__R_Wilson]

Updates:
I have better folder names.
I inserted the code to crack the feature class into multiple classes each with one feature to prevent issues overlaps.
I added FeatureToRaster to pre-rasterize the vectors before Zonal Statistics.
I added CalculateStatistics because that was the error it was throwing on the temp rasters that started this whole thing... (though arcpy.env.rasterStatistics = 'STATISTICS' might do the same thing, environmentally instead? thoughts welcome on that)

When I run it from PyScripter, hardcoded, it works wonderfully. Takes 10 minutes.

When I parameter it and load it into a Script Tool, it chokes. And it chokes on the Feature to Raster. It'll sit for over 30 minutes, right at FeatureToRaster, and not budge.

# Import modules
import arcpy

# Get parameters
indir = arcpy.GetParameterAsText(0)
zonefile = arcpy.GetParameterAsText(1)
outdir = arcpy.GetParameterAsText(5)
outfile = arcpy.GetParameterAsText(6)

# Set global environments and make scratch
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True
scratchdir = arcpy.CreateFolder_management("C:\\TEMP","ZoneScratch")
scratchdir = str(scratchdir)
arcpy.env.scratchWorkspace = scratchdir
scratchgdb = arcpy.env.scratchGDB
scratchgdb = str(scratchgdb)

# Check the input data for the presence of the zone field in zone file
# and that rasters are present in the raster folder
fieldlist = arcpy.ListFields(zonefile,"*","ALL")
fieldnames = []
for field in fieldlist:
    fieldnames.append(field.name)
if not "UID" in fieldnames:
    arcpy.AddWarning(zonefield + " field does not exist in " + zonefile)
arcpy.env.workspace = indir
rasters = arcpy.ListRasters("*")
if not rasters:
    arcpy.AddWarning(indir +" contains no rasters \nZonal statistics failed")

# Create individual feature classes for every feature in zonefile
# Allows zonal statistics to be computed for overlapping zones
zonedesc = arcpy.Describe(zonefile)
zonepath = zonedesc.catalogPath
arcpy.env.workspace = zonepath
cursor = arcpy.da.SearchCursor(zonefile,["UID"])
for row in cursor:
    rowuid = row[0]
    newfc ="scratchfc" + str(rowuid)
    whereclause = "\"UID\" = " + str(rowuid)
    featureclass = arcpy.FeatureClassToFeatureClass_conversion(zonefile,scratchgdb,newfc,whereclause)
del row
del cursor

# Create list of path names for these individual feature classes
# Allows access from a different workspace
arcpy.env.workspace = scratchgdb
features = arcpy.ListFeatureClasses("*")
featurepaths = []
for feature in features:
    desc = arcpy.Describe(feature)
    featurepaths.append(desc.catalogPath)

# Set variables and environment for zonal statistics and add field routines
arcpy.AddMessage("Setting variables and environment for ZoneStats/AddField")
zonefield = "UID"
rzonefield = "VALUE"
outtable = ""
arcpy.env.workspace = indir
arcpy.env.snapRaster = rasters[0]
arcpy.env.cellSize = rasters[0]
spatialRef = arcpy.Describe(rasters[0]).spatialReference
arcpy.env.outputCoordinateSystem = spatialRef

# Zonal Statistics
for feature in featurepaths:
    index = 0
    fname = feature.split("\\")[-1]
    zraster = scratchdir + "\\" + fname
    arcpy.FeatureToRaster_conversion(feature,zonefield,zraster)
    arcpy.AddMessage("Feature To Raster done, about to Calculate Statistics!")
    arcpy.CalculateStatistics_management(zraster)
    arcpy.AddMessage("Calculating Statistics Done")
    for raster in rasters:
        arcpy.AddMessage("Got into the raster in rasters loop!")
        outtable = scratchdir + "\\" + fname +"_" + str(index) + ".dbf"
        arcpy.sa.ZonalStatisticsAsTable(zraster,rzonefield,raster,outtable,"DATA","ALL")
        arcpy.AddMessage("Zonal statistics have been calculated for zone: " + fname + " and raster: " + raster)
        index = index + 1

[ChrisSnyder]

One idea: Before you run your last "Zonal Statistics" part of your code, try clearing the workspace:

arcpy.env.workspace = ""

Since you are now using full paths to reference your input features, output rasters, and output zone tables, I think this may be what is messing things up.

[BenSciance]

Rather than creating brand new feature classes to run zonal stats, you could select each polygon individually and run zonal stats with it.

Basic workflow would be:
1) create an index that will loop through each polygon in your fc, n=50 for example, that you want to calculate zonal stats in: 

i=1
for i in range(50):

2) convert your polygon fc with UID=1 to feature layer.

polyLyr=arcpy.MakeFeatureLayer_management(polygonFC, outputLyr, "\"UID\"={}".format(str(i)))

3) select your new feature layer by attribute.

THEN:
4) begin another for loop that will iterate through your raster grids.
5) calculate zonal stats within your selected polygon with each of your raster grids.
6) append the output zonal stats table however you'd like.

[N__R_Wilson]

Just an update:

I went ahead with hard-coded python run out of PyScripter due to time constraints. Another issue that came up was that many of my polygons were significantly smaller than my raster cell size. I created a script to resample my rasters to a 1m resolution and that has helped as well.

I also got some suggestions from Jon Bodamer (ESRI) at the User Conference. His suggestions were list comprehensions and function definitions.

When (if!) I have the time to refine this script and get it into a script tool I'll definitely incorporate the last two suggestions as well. This was definitely a good crash course in scripting and problem solving. Thanks folks!