Request for feedback on python script to increase efficiency with cursors

Question

I’m new to python and would love some suggestions on improving the efficiency of my script.

Background: My script is intended to create a report of the dominant and limiting soil types for agricultural fields. The dominant soil is simply the soil type with the greatest area in a field. The limiting soil is the soil type with greatest area, but whose name appears on the limiting soil list. So, in some fields, the limiting soil may be the same as the dominant soil, while other fields may have no limiting soils at all.

The script needs to work for all the parcels owned by an individual landowner. Each landowner may have multiple properties (called tracts—each tract has a unique ID number ), and each tract has multiple fields (numbered 1, 2, 3,… etc). The parcel data for each landowner is stored in a single feature class. It is important to be able to sort the final output by tract and field

Problem: The current script is slow! I was able to figure out how to accomplish my intended goal, but the script runs very slowly. I think it is because of my reliance on cursors and multiple if then loops. Is there a more efficient or elegant way to do this? It takes 15 to 20 min to work through some of our mid-sized landowners and I haven’t dared to try it on our largest ones.

I am running Arc 10.0 and Pythonwin 2.6.5, so I cannot use DA cursors. The main area of trouble is in the portion of code that builds the dictionaries (following "# loop through by tract"). I timed the script on a mid-sized (for this application) dataset that had 682 rows (9 tracts) and it took 00:14:16 to complete, 00:13:02 of which was in this portion of code.

import arcpy
import collections

# set environments
arcpy.env.overwriteOutput = True
workspace = r"H:\pathway..."
arcpy.env.workspace = workspace

# input data
soil = #feature class of soil polygons
parcels = #feature class of land units
outputTableMaster = #emptyDbaseTable

#intersect soil polygons and parcel layer
inFeatures = [soil, parcels]
soilTable = #output featureclass from intersect
arcpy.Intersect_analysis(inFeatures, soilTable)

# create list of limiting soils
limitingSoils = ['Au', 'Bg', 'Br', 'Ca', 'Cb', 'Ce', 'Co', 'Ru']

#identify whether or not soil types found on property are limiting
arcpy.AddField_management(soilTable, "Limiting", "TEXT", "", "", 1)
urows = arcpy.UpdateCursor(soilTable)
for row in urows:
    soilType = row.getValue("MUSYM")
    soilShort = str(soilType)[0:2]
    if soilShort in limitingSoils:
        ##soil is limiting
        row.setValue("Limiting", "Y")
    else:
        ##soil is not limiting
        row.setValue("Limiting", "N") 
    urows.updateRow(row)

del urows
del row

# create list of tracts
tractList = []
srows = arcpy.SearchCursor(soilTable)
for row in srows:
    tract = row.getValue("TRACTNBR")
    if tract not in tractList:
        tractList.append(tract)

del srows
del row

# loop through by tract and field to build dictionary
for tract in tractList:
    dictionary = collections.defaultdict(list)
    fieldList = []
    srows = arcpy.SearchCursor(soilTable)
    for row in srows:
        field = row.getValue("FieldNBR")
        if field not in fieldList:
            fieldList.append(field)
    del srows
    del row

    for field in fieldList:
        areaList = []
        srows = arcpy.SearchCursor(soilTable)
        for row in srows:
            if field == row.getValue("FieldNBR"):
                area = row.getValue("Shape_Area")
                soil = row.getValue("SoilName")
                limiting = row.getValue("Limiting")
                if area not in areaList:
                    dictionary[field].append([area, soil, limiting])
    del srows
    del row

#use dictionary to generate report of dominant and limiting soils for each field
    for data, dataList in dictionary.iteritems():
        field = str(data)
        dominantSoil = max(dataList)[1]
        domSoilLimitStatus = max(dataList)[2]
        domSoilArea = str(max(dataList)[0])
        if domSoilLimitStatus == "Y":
            limitingSoilName = dominantSoil
            limitingSoilArea = domSoilArea
        elif domSoilLimitStatus == "N":
            limitedList = [x for x in dataList if x[2] == "Y"]
            if not limitedList:
                limitingSoilName = "None"
                limitingSoilArea = "n/a"
            if limitedList:
                limitingSoilName = max(limitedList)[1]
                limitingSoilArea = str(max(limitedList)[0])

        #update the outputTable for each field         
        irows = arcpy.InsertCursor(outputTable)
        row = irows.newRow()
        row.setValue("Tract", str(tract))
        row.setValue("Field", field)
        row.setValue("DominantSoil", dominantSoil)
        row.setValue("DomSoilArea", domSoilArea)
        row.setValue("LimitingSoil", limitingSoilName)
        row.setValue("LimitingSoilArea", limitingSoilArea)
        irows.insertRow(row)

        del irows
        del row

Answer 1

I thought I would take some of my comments and add another answer. I haven't tested the answer, but I think it should work out alright. The only thing I don't know is if arcpy will remove selected features after running the initial CalculateField or not.

import arcpy
import collections


# set environments
arcpy.env.overwriteOutput = True
workspace = r"H:\pathway..."
arcpy.env.workspace = workspace

soil = #feature class of soil polygons
parcels = #feature class of land units
outputTableMaster = #emptyDbaseTable

#intersect soil polygons and parcel layer
inFeatures = [soil, parcels]
soilTable = #output featureclass from intersect
arcpy.Intersect_analysis(inFeatures, soilTable)

# create list of limiting soils
#limitingSoils = ['Au', 'Bg', 'Br', 'Ca', 'Cb', 'Ce', 'Co', 'Ru']

#here I am going to just create an expression using SQL and LIKE with the % wildcard
expression = ' "MUSYM" LIKE \'Au%\'  or "MUSYM" LIKE \'Bg%\'  or "MUSYM" LIKE \'Br%\'  or "MUSYM" LIKE \'Ca%\'  or "MUSYM" LIKE \'Cb%\'  or "MUSYM" LIKE \'Ce%\'  or "MUSYM" LIKE \'Co%\'  or "MUSYM" LIKE \'Ru%\'  '


#identify whether or not soil types found on property are limiting
arcpy.AddField_management(soilTable, "Limiting", "TEXT", "", "", 1)

#select ones to assign Y to
arcpy.SelectLayerByAttribute_management(soilTable, "NEW_SELECTION", expression)

#calculate field
arcpy.CalculateField_management(soilTable, "Limiting", "Y")

#switch selection and run again
arcpy.SelectLayerByAttribute_management(soilTable, "SWITCH_SELECTION")

#calculate field again
arcpy.CaculateField_management(soilTable, "Limiting", "N")

I have no idea if the SQL query will be faster/slower/the same as the other ways, but I am adding another option to the mix. I suppose you could set up a bigger SQL query and run one operation. I am not sure how to translate an UPDATE query into something ArcMap would be able to use though.

If you can post some timings from your experimenting, that would be pretty great for people with similar issues

Answer 2

I think a better approach in this case would be to use the calculate field method rather than building an update cursor. After you add the field, you do the following:

expression = 'is_suitable(!MUSYM!)'
codeblock = '''def is_suitable(soil_type):
    suitable_types = ['Au', 'Bg', 'Br', 'Ca', 'Cb', 'Ce', 'Co', 'Ru']

    if soil_type in suitable_types:
        return 'Y'
    return 'N'
    '''
arcpy.CalculateField_management(inTable, fieldName, expression, "PYTHON", codeblock)

That way you avoid the UpdateCursor yourself in different places, assuming you can't do the update in a single run.

As for the tractList. It would be better to use a set() in this case rather than a list. The collision lookup is O(1) then, rather than O(n). The way to do it would be:

srows = arcpy.SearchCursor(soilTable)
tracts = {row.getValue('TRACTNBR') for row in srows}
# Set comprehension - similar to the loop you created before

Use the same technique for the fields as well so you can have the two items to build your dictionary from.