ArcPy Iterative Select by Location

Question

I have a very large dataset with line features that represent a utility, and lines that come off of those utility lines showing connections to buildings. (example shown below - utility line in black and offshoots in purple.)

I developed some iterative code that takes my selected black lines, selects the offshoots attached to it, and checks some of the attributes. If the attribute checks fail, it takes the id code, appends it to a list, then later uses that list to select the purple lines and export them to their own feature. The code works fine, this is not a troubleshooting post. My questions is about efficiency. Has anyone tried something like this in the past, and is there a better way to do it? I ask because this code took about 10 minutes to run on a given study area, and we have >100,000 study areas. If you multiply that out, it would take ~19 years to run on the whole data set, which essentially renders this useless. The study area had ~130 black lines, and ~450 purple lines in it. The tool took 604 seconds to run, which is a little over a second per feature.

My thought is that running nested for loops that each call upon a Arcpy SearchCursor is very slow, and I was wondering if there way any way to get a similar effect with better performance. Are there any links to documentation or other posts? I assume this is a rather niche use case.

with arcpy.da.SearchCursor("Feature1", ['SHAPE@','ID','Attribute1','Attribute2']) as cursor1:
for row1 in cursor1:
    f_shape = row1[0]
    ID1 = row1[1]
    a1 = row1[2]
    a2 = row1[3]
    arcpy.SelectLayerByLocation_management("Feature2","BOUNDARY_TOUCHES",f_shape)
    with arcpy.da.SearchCursor("Feature2", ['SUBTYPE','ID_1','Attribute3','Attribute4','ID_2']) as cursor2:
        for row2 in cursor2:
            sub_type = row2[0]
            if sub_type <> 1:
                if a1 == 'criteria':
                    if (a2 == 'something1' or a2 == 'something2'):
                        if row2[2] <> 1:
                            found_list.append(row2[4])
                    elif a2 == 'something 3':
                        if row2[3] <> 1:
                            found_list.append(row2[4])
                if row2[1] <> ID1:
                    if row2[4] not in found_list:
                        found_list.append(row2[4])

Answer 1

Thanks to the suggestion from @ziggy, I realized you can preform geometry operations on features outside of an arcpy.da.Cursor by storing the '@SHAPE' attribute in python. You then simple pull in a geometric "Method" - in this case shape1['@SHAPE'].touches(shape2['@SHAPE']).

So the code ended up looking like this, which compared to my original, is drastically different. Basically, I call in a search cursor once for each feature class, storing all of the relevant values in a dictionary (which itself was an entry in another dictionary based on unique id codes - a dictionary of dictionaries).

I then did all of the logical and spacial checks directly on those dictionaries and write out the unique ids of any features that failed to a list. Eventually that list is converted to a tuple to for a select by attribute function to export out only the features that are flagged.

#create list to hold found lines - populate with 2 dummies for tuple syntax
found_list = ['dummy0','dummy1']
#create dictionary to store search cursor data
f1_dict = {}
with arcpy.da.SearchCursor(feature1, ['SHAPE@','GLOBAL','att11','att12']) as cursor1:
    for row1 in cursor1:
        f1_dict[row1[1]] = {'f_shape' : row1[0],
                              'global' : row1[1],
                              'att_11' : row1[2],
                              'att_12' : row1[3]
                              }
#create dictionary to store search cursor data
f2_dict = {}
with arcpy.da.SearchCursor(feature2, ['SHAPE@','SUB_TYPE','att21','att22','att23','GLOBAL']) as cursor2:
    for row2 in cursor2:
        f2_dict[row2[5]] = {'f_shape' : row2[0],
                              's_type' : row2[1],
                              'm_id' : row2[2],
                              'att_22' : row2[3],
                              'att_23' : row2[4],
                              'global' : row2[5]
                              }
#spatial check then write out to list
for x1 in f1_dict:
    a1 = f1_dict[x1]
    b1 = a1['f_shape']
    for x2 in f2_dict:
        if f2_dict[x2]['s_type'] <> 1:
            a2 = f2_dict[x2]
            b2 = a2['f_shape']
            if b2.touches(b1) == True:    #<<<the select by location equivalent
                if a1['att_11'] == 'ABV60':
                    if (a1['att_12'] == 'LTRAN' or a1['att_12'] == 'RTRAN'):
                        if a2['att_22'] <> 1:
                            found_list.append(x2)
                    elif a1['att_12'] == 'DIST':
                        if a2['att_23'] <> 1:
                            found_list.append(x2)
                if a2['m_id'] <> a1['global']:
                    if x2 not in found_list:
                        found_list.append(x2)


arcpy.SelectLayerByAttribute_management(feature1, "CLEAR_SELECTION")
arcpy.SelectLayerByAttribute_management(feature2, "CLEAR_SELECTION")

arcpy.SelectLayerByAttribute_management(feature2, "NEW_SELECTION", '"GLOBAL" IN '+str(tuple(found_list)))
arcpy.CopyFeatures_management(feature2, gdb1+'\\'+output1)

The end result was that from a 604 second run-time for the original code, it dropped to 36 seconds in this code, a 96% decrease. I also liked the answer provided by @Emil Brundage but a spatial join creates an intermediary file, that then has to be read in and operated on. All in all I found this solution to be the cleanest.

Answer 2

Make use of a spatial join instead of selecting each feature one at a time. This will give you a table with your main fields and your service fields based on the indicated spatial relationship. This eliminates the need for multiple spatial selections as well as multiple cursor initiations, both of which take a lot of time.

arcpy.SpatialJoin_analysis ("Feature2", "Feature1", r"some/out/featureclass", 
                            match_option = "BOUNDARY_TOUCHES")