[GIS] faster way to add a field based on zonal statistics to an attribute table

arcgis-desktoparcmapzonal statistics

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.

Best Answer

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
Related Question