Iterating Intersect_analysis over polygon shapefile using ArcPy?

Question

I am using ArcGIS 10.3.1 for Desktop and I have two shapefiles:

1. A set of stream (vector line) data
2. A fishnet of polygons

I want to know, for each polygon in that fishnet, what the total length of streams is within that polygon. So far, I've tried iterating the intersect_analysis tool over the fishnet table, summing the lengths of the output and writing that to the relevant field, but each iteration performs the intersect on the entire grid, not the selected polygon in that grid.

Code:

import arcpy
#set some environmental variables
arcpy.env.workspace = "C:\Users\Alex\Documents\SoilMoisture"
arcpy.env.overwriteOutput = 1
#define variables for all necessary feature classes and fields
streams = "pathto/str_split_lamb.shp" ## shapefile of stream network split at vertices, in Lambert projection, with Length value field pre-calculated
rainGrid = "pathto/net_test.shp" ## polygons for rain cells, with empty "DRN_LENGTH" field
cursor = arcpy.UpdateCursor(rainGrid) ## will iterate across every row of the rainfall grid feature
for row in cursor:
 intStreams = arcpy.Intersect_analysis([streams,rainGrid],"output.shp","ALL","","LINE") ## intStreams becomes all stream segments inside polygon
 strCursor = arcpy.SearchCursor(intStreams) ##iterate over streams in intersection
 totalLen = 0
 for rowmore in strCursor:
 totalLen += rowmore.length ## sum all length values into one variable
 print totalLen # this is just here to check what's going on
 row.DRN_LENGTH = totalLen ## add value to field
 cursor.updateRow(row)
print "Done!"

Any ideas?

• FelixIP: Why on earth to iterate through each polygon if you can intersect all in 1 go, summarise length by polygon id and bring this numbers back to your fishnet dataset gis.stackexchange.com/questions/169608/…

Answer 1

You should not iterate any SQL selections at all, since that is very slow and totally unnecessary. Instead intersect all lines by all polygons, then use a search cursor on the intersect results to build a dictionary with the polygon FID as the dictionary key with values that contain the sums of the line lengths. Finally use an update cursor to write to each polygon after looking up the polygon FID dictionary key. Once the Intersect is complete, the actual process of summarizing the intersect data for the full grid and writing the lengths to the target dataset will finish in seconds for 10,000 records.

See my Blog on Turbo Charging Data Manipulation with Python Cursors and Dictionaries. In particular look at the section on "Using a Python Dictionary Built using a da SearchCursor to Replace a Summary Statistics Output Table" at the end of the Blog.

• Alex Wheatley: Interesting. I wasn't aware that this would work. As it happened, I had a weekend to let it run, so in the end adding 2 lines to the code and letting it think was quicker for me, but I'll definitely try this next time I need to do something similar.

Answer 2

You need to select a single polygon with each iteration first; then the intersection tool will only consider that one polygon. You can select by attribute (http://pro.arcgis.com/en/pro-app/tool-reference/data-management/select-layer-by-attribute.htm), using unique FIDs.

• Alex Wheatley: Aha...I hadn't quite understood exactly what the cursor was doing. That makes a lot of sense.

Answer 3

Building on @Paulo Raposo answer, if the script will be run outside of the Python window in ArcMap and the layer is not in the current map, you will need to use Make Feature Layer method before performing feature selection using Select Layer by Attribute.