Selecting all polygons where neighbor vertex points are below distance threshold in ArcPy?

Question

I have a shapefile with over 100k polygons (some of them have holes in them but apart from that they are simple polygons - no multipart). I want to loop through all of these polygons and select only the ones where at least two neighbor vertices (within the same polygon) are within a distance threshold. Another way of putting it is: if two neighbor vertices (ex: id #3 and id #4) within a polygon have a distance < threshold_distance, then select them for further processing. Any idea how to do this efficiently?

I am a beginner with arcpy (python). Conceptually, the brute force approach I see is:

1. Use a search cursor to loop through all polygons
2. For each polygon, use "getPart" to create an array of points
3. Use a for loop to cycle through all points and calculate distance between point i and point i+1 (except for last point). If the distance is less than distance_threshold, then stop the loop, copy the polygon and insert it into a new shp file. (or maybe just copy the polygon's unique identifier and save that into a txt file)

Since I need to cycle through + 100k polygons, is there a more efficient way of doing things?

Answer 1

Efficiency might not be my strong suit, but here is the process I came up with:

In the code below, I first create a point feature class from the vertices of my input shapefile. From here I iterate through each feature in the input shapefile. I select each vertex in the point feature class associated with that feature, and then perform a distance selection. After unselecting the original feature points, I iterate through my remaining point selection and check if two are selected with the same original FID field values. If so, I have indicated where in the code you can perform further geoprocessing. I'm curious what further geoprocessing you are planning on implementing. I might be more efficient to do so outside of the nested loop.

Below, 'inputFC' is a string of the full path of your shapefile. Make sure you have the 'r' before your open quotes. 'distance' is the distance you are using.

Code:

import arcpy

inputFC = r"C:\GISStuff\input.shp"
distance = "100 Meters"

#Create points from input feature
arcpy.FeatureVerticesToPoints_management (inputFC, r"in_memory\points")

#Create layer files to allow selection
arcpy.MakeFeatureLayer_management (r"in_memory\points", "pntlyr")

#Create cursor for table access
cursor = arcpy.da.SearchCursor (inputFC, ["OID@"])
for row in cursor:
    #list to append first point per OID to
    onePntli = []
    #list to append second point per OID to
    twoPntli = []
    #SQL to select points by original OID
    sql = '"ORIG_FID" = ' + str (row[0])
    #Select points for feature
    arcpy.SelectLayerByAttribute_management ("pntlyr", "", sql)
    #Select neighboring points
    arcpy.SelectLayerByLocation_management ("pntlyr", "WITHIN_A_DISTANCE", "pntlyr", distance)
    #Remove original feature points from selection
    arcpy.SelectLayerByAttribute_management ("pntlyr", "REMOVE_FROM_SELECTION", sql)
    #Check for selection:
    if arcpy.Describe ("pntlyr").FIDSet:
        #cursor to iterate through selected points
        pntcursor = arcpy.da.SearchCursor ("pntlyr", ["ORIG_FID"])
        for pntrow in pntcursor:
            #If one point has been found but not two...
            if pntrow[0] in onePntli and not pntrow[0] in twoPntli:
                #Add ORIG_FID to list indicating two points have been found
                twoPntli.append (pntrow[0])
                ###
                #
                ###!!!Further Geoprocessing here!!!
                ### variable pntrow[0] contains the OID of the neighboring polygon  
                ### feature to be geoprocessed
                #
                ###
            #If this is the first point found with given "ORIG_FID"...
            else:
                #Add pnt "ORIG_FID" to list
                onePntli.append (pntrow[0])
        #Clean up
        del pntrow
        del pntcursor
del row
del cursor
arcpy.Delete_management ("pntlyr")
arcpy.Delete_management (r"in_memory\points")

I hope this helps.