1

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.

  • What is wrong with that workflow? It seems fine to me. Have you tried just 'Zonal Statistics' and not Zonal Statistics as Table? From memory this adds a field to the polygon table with the statistic type but you have no control over what the field is called. – Michael Stimson Feb 15 '17 at 23:45
  • @MichaelMiles-Stimson Mainly it just cumbersome to do repeatedly, joining and unjoining and generating a table I don't need for anything else. Unless I am mistaken, zonal statistics just produces a raster of the zone values. – K. Credo Feb 15 '17 at 23:58
  • I'm imagining something like Extract Values to Points, except this would extract zonal averages to features. – K. Credo Feb 15 '17 at 23:59
  • There are now other single click tool. Create a model, convert to script. Replace joining/field calculator by dictionary/da.updatecursor and will work like a charm. Alternatively reuse lines in Results window – FelixIP Feb 16 '17 at 3:25
  • That's the way I'd do it @FelixIP, digest up to a dictionary and dish out with an update cursor - works well and less likely to crash out randomly. DBF files must be less than 2GiB so it's unlikely to cause memory issues in a separate (32bit in 64bit) memory space but is something to consider if you're on a dodgy notebook with pitiful RAM and lots of polygons. – Michael Stimson Feb 16 '17 at 4:55
3

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
  • 1
    +1 nice. I'd shorten statistics table name, just in case. And da cursor is fundamentally faster than field calculator – FelixIP Feb 16 '17 at 5:04
  • @FelixIP have we ever performance tested to validate your statement that "da cursor is fundamentally faster than field calculator"? I would have expected the opposite, but with there not being a huge difference, unlike between old and new style cursors. – PolyGeo Feb 16 '17 at 5:14
  • 1
    @PolyGeo according to gis.stackexchange.com/questions/198163/… calculator twice slower. 200k were used. At 1000k the difference even greater. I cannot find my own answer, where surprise surprise arcview 3 beat them all – FelixIP Feb 16 '17 at 7:04
  • @FelixIP Sometime I should test what I would call fundamental field calculations (field = 3, field1 = field2) on attribute table rather than joined fields. I'm not sure but I think they would stay in C and away from the Python Parser. – PolyGeo Feb 16 '17 at 7:35
1

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.

  • +1 nice and easy to read, Please fix tab for except block – FelixIP Feb 16 '17 at 22:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.