Building A Soils Clip Tool for Agricultural Preservation Reporting

I’m getting close to wrapping up a project for our Planning department and wanted to share about it.  It’s a custom ArcGIS Desktop geoprocessing tool that stemmed from a conversation with one of our planners.  He was having trouble clipping a dataset.  After figuring out the issue, I had him walk me through his entire workflow.  Throughout our discussion, it became clear that there were several manual tasks that could be automated by a custom tool.

As part of our County’s agricultural land preservation program, incoming applications (farms) are processed and scored.  This process includes clipping the soils dataset to the farm, and calculating the total acres of soil by both soil name and soil class.  While my initial goal was to develop one tool that could handle two approaches to the input farm layer, I developed two tools to accomplish this.

Soils clip tool that takes parcel Id’s as the farm input

In most cases, the user will enter one or more parcel identification numbers that represent the farm. But occasionally, they will create a layer representing the farm.  This is because a portion of the parcel will not be included in the agricultural easement.  I learned that having a single tool that uses conditional logic to allow the user to select between the parcels or farm layer inputs would be impossible after posting on Geonet.

While the farm input is different between the two tools, everything else is the same.  Here are the user inputs for the tool:

1. Name of the farm
2. Surveyed or deeded acres of the farm
3. Parcel ID’s or layer representing the farm
4. Directory to place output Excel files creating by the tool

Soils clip tool that takes a GIS layer as the farm input

The first step in the tool is to clip a hard-coded soils dataset to the farm.  For the farm layer tool, I pass in the farm GIS layer as the layer to clip the soils against.  But for the tool the constructs the farm based upon parcel IDs, there’s a litttle more work.

I loop through the list of parcel ID’s and perform a series of select by attribute operations.  To make sure each selection actually selects a record, I perform a test where I compare the number of selected records from the previous loop iteration to the current loop iteration.  If the number is greater, I keep iterating. If not, I add a warning message to the user and exit the script.  While this helps ensure all parcel ID’s entered check against the input parcels layer, it can’t stop the user from mis-entering a parcel ID that exists in the dataset (false positive).

# make feature layer for parcels
# taxParcels is a variable representing the tax parcels dataset
arcpy.MakeFeatureLayer_management(taxParcels,'selectedParcels')

# counter variable for parcels list to test if more than one feature is selected if multiple parcels are entered
selectionCount = 1

# pinList is obtained from the tool's user input form
# parcel represents a parcel ID number
for parcel in pinList.split(';'):
        # where clause
        whereClause = "APPENDPIN = '{}'".format(parcel)
        # creates a new selection for first parcel
        # adds to existing selection for all others
        arcpy.SelectLayerByAttribute_management('selectedParcels', 'ADD_TO_SELECTION', whereClause)
        # get count of parcels selected in loop
# convert to integer to use in conditional test
        matchCountParcels = int(arcpy.GetCount_management('selectedParcels')[0])
        # test if number of selected parcels has changed since previous iteration; skip during first selection
        if selectionCount > 1:
            if matchCountParcels > numberSelectedParcels:
                # do not do anything if more parcels are selected
                pass
            # if number of selected parcels has not changed
            else:
                # add a warning that the entered parcel number did not result in additional selections
                arcpy.AddWarning('\nNo parcel(s) were selected for APPENDPIN #{}.  Please verify the APPENDPIN numbers entered in the form and run the tool again\n'.format(parcel))
                # exit tool
                sys.exit()
        # number of selected parcels from current iteration
        numberSelectedParcels = matchCountParcels
        # increase counter
        selectionCount += 1
    # end for

Now, let me outline how the rest of the tool flows:

1. Soils are clipped to farm. Output layer is given standardized name and placed in a hard-coded geodatabase.
2. A field is added to represent the ratio acres for each soil feature
3. The total acres of the clipped dataset is calculated using arcpy.da.FeatureClassToNumPyArray() and the sum() method on the NumPy array
4. A constant is created by dividing the surveyed/deeded acres by the summed acres (#3).
5. The ratio acres field created in #2 is calculated for each record by multiplying the ACRES field by the constant (#4)
6. The clipped soils dataset is exported to a Microsoft Excel file and placed in a directory specified by the user
7. A table is created in memory to store summarized data for the clipped soils dataset based upon soil name.  The table includes the soil name and soil class fields
8. A field is created to store the percent of total value for each feature
9. The total ratio acres are summed using a similar method to #3.
10. The percent of total field (#8) is calculated for each feature by dividing the ratio acres field by the sum of the ratio acres
11. The in-memory table (#7) is exported to a Microsoft Excel file and placed in a directory specified by the user
12. Steps 7-11 are repeated, but based upon the soil class field
13. The in memory workspace is deleted
14. The three output Excel files include the clipped soils dataset, total ratio acres by soil name, and total ratio acres by soil class

Preserved farms (green) in Cumberland County, Pennsylvania

The last topic I wanted to cover was an error handling function I’ve incorporated in this tool.  And honestly, I plan to use it in all of my custom tools.  I first found this function in a tool I found on Geonet.  That tool allows you to copy data from map and feature services to a local dataset.

If an error occurs during the tools execution, an error message appears that directs the user to provide information to our department.  The information includes the system error message, python file where the error occurred, and line number where the error occurred.  I’m hoping this will make it easier to diagnose and fix issues that will inevitably arise as this tool is used.  Here’s the snippet:

# This code is in a separate script I import into the tool's script
# errorLogger.py

# import modules
import sys, linecache, arcpy

# Function to handle errors
def PrintException(error):
    exc_type, exc_obj, tb = sys.exc_info()
    f = tb.tb_frame
    lineno = tb.tb_lineno
    filename = f.f_code.co_filename
    linecache.checkcache(filename)
    line = linecache.getline(filename, lineno, f.f_globals)
    arcpy.AddError('\nerror: {}\nFILE: {}, LINE: {}\n\n\t "{}": {}'.format(error, filename, lineno, line.strip(), exc_obj))
    # exit Python
    sys.exit()
# end PrintException

# And then within tool's script
# soilsClipTool.py
import errorLogger, arcpy # and the other modules
try:
   # all the tool's logic
   pass
# If an error occurs running geoprocessing tool(s) capture error and write message
# handle error outside of Python system
except EnvironmentError as e:
    arcpy.AddError('\nAn error occured running this tool. Please provide the GIS Department the following error messages:')
    # call error logger method
    errorLogger.PrintException(e)

# handle exception error
except Exception as e:
    arcpy.AddError('\nAn error occured running this tool. Please provide the GIS Department the following error messages:')
    # call error logger method
    errorLogger.PrintException(e)
# delete in_memory workspace
finally:
    try:
        # attempt to delete in_memory workspace
        arcpy.Delete_management("in_memory")
    except:
        pass