2

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?

2

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.

  • Thanks for the input! Once I have identified all polygons in which there are closely spaced vertices, I must then eliminate one of the two vertices (a kind of simplification of the polygon). Which vertex to eliminate will depend on the topology of the polygon with respect to the neighboring polygons and will have to be done on a case by case basis (eg manually!) – lavaman Dec 23 '14 at 18:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.