1

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
  • 1
    Welcome to GIS SE. Can you edit your post/tags to include the version of ArcMap that you are using? – Barbarossa Aug 6 '15 at 17:43
  • 1
    Have you considered trying a selection query (help.arcgis.com/EN/arcgisdesktop/10.0/help/index.html#//…) and then running a field calculator (help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//…) on the selected features? I am not sure how speed will be, but running a SQL query will probably be quicker than a cursor... You can also look into the da version of the cursor (arcpy.da.UpdateCursor) as it is built more for speed. – Branco Aug 6 '15 at 17:53
  • Yes, the data access (da) cursors are much faster than the older cursors. However, the da module was not introduced until ArcMap 10.1. – Barbarossa Aug 6 '15 at 18:08
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

  • @Barbarossa.. that is some pretty interesting information. Guess I probably should have searched knowing someone, somewhere had to of done a test or two. – Branco Aug 6 '15 at 19:47
  • I'm a big advocate of da cursors and I try to use them exclusively. I find huge gains over other cursors and field calculations. For instance, I rewrote an old script that used SQL queries and selections to use da cursors. Processing time in ArcMap literally went from 10 min to 10 sec. – Barbarossa Aug 6 '15 at 19:48
  • Granted, if your dataset is relatively small, either way should be sufficiently fast. However, when you move into tens of thousands of records, speed gains are nice. – Barbarossa Aug 6 '15 at 19:55
  • I just never get anything larger than maybe 10-500 features so it really doesn't matter for me. I typically stick with the regular cursors when editing data for ease of use. Plus, most of my scripts go to analysts/specialists where they look through and tweak so it is easier for them to understand. Give me a large data set with some real work to be done, I will definitely go the da route. – Branco Aug 6 '15 at 19:55
0

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.

  • thanks for the suggestion! I'm new to sets and having trouble with the syntax you provided. Does this construction work with my version of Arc/Python? Do I need to import a module to start? or could there be a typo in my code? tracts = {row.getValue("TRACTNBR") for row in srows} – gego Aug 7 '15 at 12:37
  • Maybe there is an issue with your version. I don't have arcpy installed so I can't check. The long version would be: tracts = set(); for row in srows: tracts.add(row.getValue('TRACTNBR')); – Michalis Avraam Aug 20 '15 at 23:26

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.