Using ArcPy Search Cursors to Perform Spatial Analysis by Each Record in a Jurisdictional Dataset

road closed near a river where flooding occured

Cumberland County is in the process of updating our Hazard Mitigation Plan.  Our GIS Department is focused on developing a web map series application, and performing various spatial analyses to support the plan.

One of the first analyses we were asked to complete was to get a count of the number of structures in the floodplain, by municipality. As this would require running the same steps over and over again, I thought that writing a Python script that uses a for-in loop would save some time.  And through this process, I was able to develop a more generic, re-usable geoprocessing tool. You can download the tool on GitHub.

This tool is great whenever you need to get all of the features from Dataset A (i.e., building structures) that intersect Dataset B (i.e., floodplains), for each feature in Dataset C (i.e., municipalities).  Performing this type of analysis manually could be time-consuming.  But creating a tool to automate these repeated steps saves you time now, and in the future when you complete a slightly different project.

In this article, I’ll review the Python script for this specific building structures in floodplain analysis.  I’ll also share the geoprocessing tool I developed out of this process.  Hopefully you can use the tool in one of your future projects.  Let’s begin!

road closed near a river where flooding occured
A roadway is closed along the Susquehanna River in Enola, PA.

Looping Through The Dataset Using a Search Cursor

Because we want a separate list of buildings that intersect floodplains for each municipality, we’ll utilize a for-in loop using a Search Cursor (arcpy.da.SearchCursor) on the municipalities dataset.  Below is a simplified diagram of the analyses that occur within each iteration of the loop.  I’ll also outline these steps using a list:

  1. Create an attribute selection on a feature layer version of the municipality dataset using the current record in the loop. (feature layers are required for attribute or spatial selections in ArcPy)
  2. Select features from the building structures dataset that have their centroid in the selected municipality feature
  3. Execute code based upon whether any building structures are selected or not:
    1. Features are Selected:
      1. If features are selected, select features from the building structures (subset selection) that intersect the floodplains dataset.
      2. Export selected records to a new dataset in the project folder or geodatabase using the current municipality in the file name.
      3. If no building structures intersect the floodplains dataset, add a message stating this.
    2. No Features are Selected:
      1. Add a message that no building structures have their centroids in the municipality.
diagram of python script that selects building structures that intersect floodplains for each municipality
Diagram of script. The series of selections are performed on each record in the Municipalities dataset.

Below is the sample code for the Search Cursor portion of the script:

# import arcpy and os modules
import arcpy, os
# Create text file for logging results of script
logFile = r'[path]\[to]\[folder]\Buildings_Intersect_Floodplains_Report.txt'
# variable to store messages for log file. Messages written in finally statement at end of script (not shown in sample)
logMsg = ''
# output location for shapefiles
out_folder_shp = r'[path]\[to]\[folder]\Buildings In Floodplains\Geodata'
# building footprints
bldg_footprints = r'[path]\[to]\[folder]\Geodata.gdb\Building_Footprints'
# fema floodplains
fema_floodplains = r'[path]\[to]\[folder]\Geodata.gdb\Floodplains'
# municipalities
muni = r'[path]\[to]\[folder]\Geodata.gdb\Municipality'
# municipalities fields
muni_fields = ['MUNI_NAME']
# create feature layer for municipalities
# create feature layer for buildings
arcpy.MakeFeatureLayer_management(bldg_footprints, 'buildings')
# create feature layer for floodplains
arcpy.MakeFeatureLayer_management(fema_floodplains, 'floodplains')
# Create Search Cursor for Municipalities
# add message
logMsg += 'Creating search cursor on {}\n'.format(muni)
# create search cursor
with arcpy.da.SearchCursor(muni,muni_fields) as cursor:
# loop through each record in municipalities layer
for row in cursor:
# add message
logMsg += '\nSelecting buildings that intersect {}\n'.format(row[0])
# attribute selection on municipality
# SQL query
where_clause = "MUNI_NAME = '{}'".format(row[0])
arcpy.SelectLayerByAttribute_management('municipality', 'NEW_SELECTION', where_clause)
# select buildings that have centroid in municipality
arcpy.SelectLayerByLocation_management('buildings', 'HAVE_THEIR_CENTER_IN', 'municipality', selection_type='NEW_SELECTION')
# get count of buildings that intersect municipality
match_count_bldg = int(arcpy.GetCount_management('buildings')[0])
# make sure features are selected before moving forward
if match_count_bldg > 0:
# select from selected buildings that intersect floodplains
arcpy.SelectLayerByLocation_management('buildings', 'INTERSECT', 'floodplains', selection_type='SUBSET_SELECTION')
# get count of selected buildings
match_count_blg_floodplains = int(arcpy.GetCount_management('buildings')[0])
# make sure features are selected before moving forward
if match_count_blg_floodplains > 0:
# add message – number of buildings in municipality
logMsg += '\nThere are {} buildings that intersect floodplains in {}\n'.format(match_count_blg_floodplains, row[0])
# export layer
arcpy.CopyFeatures_management('buildings', os.path.join(out_folder_shp, 'Buildings_Intersect_Floodplains_{}.shp'.format(row[0])))
logMsg += '\nNo buildings that are within {} intersect floodplains\n'.format(row[0])
logMsg += '\nThere were no buildings that have their centroids in {}\n'.format(row[0])
# end if/else
# end for
# end with cursor

Building the Geoprocessing Tool

Though I wrote this Python script for a specific purpose, I soon realized I could modify it to create a re-usable tool. This tool can be used whenever you need to get all of the features from Dataset A that intersect (or another spatial selection type) Dataset B, for each feature in Dataset C.

You can download the tool on GitHub.  There could definitely be some improvements to the messages created for users.  And maybe someday, I’ll make those (or maybe you can and submit a pull request!).  I have some screenshots of the tool below.  I’ve also listed a sample use case for the tool.  I’m hoping you can use/modify/improve the tool for use with your organization.

screen-shot of a custom ArcGIS geoprocessing tool's user form

The user form for the custom geoprocessing tool

Sample Use Case

Select all Class A Wild Trout Streams that are within [x] distance of an Oil/Gas well, by County.  When performing a spatial selection that involves a distance, it is important to remember that the “tool evaluates a spatial relationship in the coordinate system of the Input Feature Layer data source (the feature class on disk)” (see Esri help).  It is possible to set the output coordinate system to evaluate the spatial relationship in.  That is another potential upgrade to this tool.

Here are the links to the datasets for this use case:

screen-shot of system messages for an ArcGIS custom geoprocessing tool

Messages from the custom geoprocessing tool

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s