Is there a faster way to add a field based on zonal statistics to an attribute table?

Question

Let's say I have a shapefile 'Mypolygons' with an ID field and a raster 'DEM' giving elevation values. I want to add a field for each row in the 'Mypolygons' attribute table to give the average elevation within each polygon.

My current approach is the following:

• Add field "Avg_Elev" to 'Mypolygons'.

• Use Zonal Statistics as Table with 'Mypolygons' as the zone data, 'ID' as the Zone field and 'DEM' as the in_value_raster, and 'Output' as the out_table to hold mean elevation values.

• Join 'Mypolygons' and 'Output' based on the 'ID' field.
• Use Calculate field to populate values of 'Mypolygons' 'Avg_Elev' field based on the 'MEAN' field in 'Output'.
• Remove join.

Is there a single tool or a more efficient way to perform this task? A python solution will work for me as well.

Answer 1

Your workflow is adequate and would easily be incorporated into a model but for a few tricky spots with joins...

Qualified field names are a pain in the join 'to' table, you can turn them off but they will suddenly fully qualify if the field name already exists in the from table (usually at a minute to knock-off on Friday). To avoid this inconsistency take charge of the temporary table name so you can predict the fully qualified field, in this case '!RasterStatTable.{}!'.format(StatType) - I have made the stat type a variable so as to make the code reusable for other statistic types but beware some fields in the table are not named the same as their statistic type.

Or alternately not use a join at all as @FelixIP said in his comment, this is the way that I'd go - not any quicker but much less likely to fail randomly. Digest all the table into a dictionary object (you should have enough RAM to load the entire table) then update the polygons using the dictionary.

I've put both methods into this code, try it for yourself and see what you prefer:

import os, sys, arcpy

PolygonFC = sys.argv[1]
Raster    = sys.argv[2]
StatType  = "MEAN"
ToField   = "Avg_Elev"
Temp      = os.environ.get("Temp")
TempTab   = os.path.join(Temp,"RasterStatTable.dbf")

# get the OID field from the feature class, that way we can work with
# shapefiles and geodatabase feature classes
Desc      = arcpy.Describe(PolygonFC)
OID_Field = Desc.OIDFieldName         # FID for shapefile, OBJECTID for geodatabase

# clear the path for the temp table
if arcpy.Exists(TempTab):
    try:
        arcpy.Delete_management(TempTab)
    except:
        arcpy.AddError("Unable to clear temp table")
        sys.exit(-1)

# get spatial analyst extension, required for the tool
if arcpy.CheckExtension("Spatial") == "Available":
    arcpy.AddMessage( "Checking out Spatial")
    arcpy.CheckOutExtension("Spatial")
else:
    arcpy.AddError("Unable to get spatial analyst extension")
    sys.exit(-2)

arcpy.AddMessage("Performing Zonal Statistics as Table")
arcpy.sa.ZonalStatisticsAsTable(PolygonFC,OID_Field,Raster,TempTab,'DATA','MEAN')

# check if the 'to' field exists and if not found then add it
if len(arcpy.ListFields(PolygonFC,ToField)) == 0:
    arcpy.AddMessage("Adding Avg_Elev field")
    arcpy.AddField_management(PolygonFC,ToField,"DOUBLE")

CalcByJoin = False # change this switch to update the table by join/calculate or dict/update
if CalcByJoin:
    # joins are annoying but you #should# be able to do it this way
    # joins must be performed on a Layer or Table View object...
    arcpy.MakeFeatureLayer_management(PolygonFC,"Layer")
    arcpy.AddJoin_management("Layer",OID_Field,TempTab,"FID_")
    arcpy.CalculateField_management("Layer",ToField,'!RasterStatTable.{}!'.format(StatType),"PYTHON")
    arcpy.RemoveJoin_management("Layer","RasterStatTable") # not really necessary, it dissapears at the end of the script anyway
else:
    # this works too and is a lot less prone to random errors
    # digest the table into a dictionary
    arcpy.AddMessage("Digesting table")
    ZonStatDict = {}
    with arcpy.da.SearchCursor(TempTab,["FID_",StatType]) as SCur:
        for ThisRow in SCur:
            ZonStatDict[ThisRow[0]]=ThisRow[1]

    # Update the feature class with the digested table dictionary
    arcpy.AddMessage("calculating values")
    with arcpy.da.UpdateCursor(PolygonFC,[OID_Field,ToField]) as UCur:
        for ThisRow in UCur:
            if ThisRow[0] in ZonStatDict:
                ThisRow[1] = ZonStatDict[ThisRow[0]]
                UCur.updateRow(ThisRow)
            else:
                # every key (OID) should exist in the dict but in SDE anything can happen in between...
                arcpy.AddWarning("Key {} not found, no value calculated".format(ThisRow[0]))

arcpy.Delete_management(TempTab)  # remove temp table, just for neatness.. if the join still exists this wouldn't work.
arcpy.CheckInExtension("Spatial") # return spatial analyst to the pool

Answer 2

Many thanks for all the great information and suggestions. It is useful and valuable to know from more knowledgeable folks that there is not an easier way to perform this task.

Here's the simplified Python version I came up with, based on code by @MichaelMiles-Stimson and all your comments:

features = "path_to_features"
zone_field ="ID_field"
input_raster = "My_DEM"
field_name = "Avg_Elev"
field_type = "DOUBLE"
stat = "MEAN"

import arcpy
from arcpy.sa import *

def add_zonal_field(features, zone_field, input_raster, field_name, field_type, stat):

    """ Performs zonal statistics on a set of features, 
    then adds and populates a new field in the feature class 
    to store the results of the zonal calculations"""

    # Add zone field to features
    arcpy.AddField_management(features, field_name, field_type)

    # Clear path for temporary table
    if arcpy.Exists("zonal_table"):
        try:
                arcpy.Delete_management("zonal_table")
        except:
                arcpy.AddError("Unable to clear temp table")
                sys.exit(-1)


    # Zonal statistics
    arcpy.AddMessage("Performing zonal statistics: " + field_name + " in " + features)

    zonal_table = ZonalStatisticsAsTable(features,zone_field, input_raster, "zonal_table", "DATA", stat)


    # Digest statistics from zonal_table
    arcpy.AddMessage("Digesting " + stat + " from zonal table")

    stat_dict = {}

    with arcpy.da.SearchCursor(zonal_table,[zone_field,stat]) as cursor:
            for row in cursor:
            stat_dict[row[0]] = row[1]


    # update new field in feature class
    arcpy.AddMessage("Calculating " + field_name + "in " + features)

    with arcpy.da.UpdateCursor(features, [zone_field, field_name]) as cursor2:
        for row2 in cursor2:
            row2[1] = stat_dict[row2[0]]
            cursor2.updateRow(row2)

I have not heard of this term "digesting" before, but to me this is very preferable to using Join. It may be my amateur approaches, but I seem to always run into problems trying to automate Joins, whether in Python or ModelBuilder. I greatly appreciate the feedback on this, as it is extremely useful for progressing from rudimentary coding to more polished versions.