region growing feature selection using arcpy

Question

I have a polygon featureclass consisting of parcels touching boundary with each other and for which areas are computed in the attribute table (Figure a). I want to select polygons in an iterative way until a given threshold of area is reached.

In my system, I want to select any one of the parcels (For example, the smallest one) as a first step (Figure b). As a second step I want to apply region growing approach that will select a neighboring parcel and add it to the present selection (Figure c and d). However, when first three parcels are selected, their total area is equivalent to 1234 unit. I want my algorithm to stop when the area reaches , for example, 942 unit. For this, I want to split the latest selected parcel to equalize the threshold (figure e).

I am familiar with basic pyhton and arcpy stuff. If someone could give me the framework to implement this.

data = "C:\Users\ln88615\Documents\GIS_data\cadastral_computation\mask2.shp"
CadastreOut = "C:\\Users\\ln88615\\Documents\\GIS_data\\cadastral_computation\\"
mask_layer = "mask_layer"
arcpy.MakeFeatureLayer_management(data, mask_layer)

field = "Matr"
cursor = arcpy.SearchCursor(mask_layer)

for row in cursor:
   geometry = row.getValue("Shape")
   out = str(row.getValue(field))
   arcpy.SelectLayerByLocation_management (mask_layer, "BOUNDARY_TOUCHES",geometry)
   arcpy.CopyFeatures_management(mask_layer, CadastreOut + "Cadastre_" + out + ".shp")

Answer 1

Just completed working on a very similar problem to yours to grow a region. Given a study area full of parcels, I needed to create polygons covering the study area combining parcels into groups of at least 1000 and no more than 5000. The code is not necessarily the most efficient but it does get the job done and may be helpful to you.

In the code below, Join_Count would be analogous to your shape area. The general idea is to select the touching neighbours and then if the sum amount exceeds the threshold you start removing features to get back under the threshold.

In your situation, once you get back under the threshold you would then figure out how much you need to cut the polygon to add back the area to reach your threshold.

I should also add that the code below is completely automated with no user interaction. I simply start with the lowest count and work up.

#---------------------------------------------------------------------------------------------------------
# Need to aggregate polygons with small counts to neighbours.
# Since the feature class needs to be updated at the same time as we work through it we will continuously
# open and close cursors to get the smallest count and then delete features that get merged to it.
#---------------------------------------------------------------------------------------------------------
geom = arcpy.Geometry()
flds = ['SHAPE@', 'OID@', 'Join_Count']
wcCount = '"Join_Count" < ' + str(minThreshold)

arcpy.MakeFeatureLayer_management(fcSJ, lyrSJ)
with arcpy.da.SearchCursor(fcSJ, flds, wcCount) as cur:
    for row in sorted(cur, key=operator.itemgetter(2)):

        #------------------------------
        # Get the features that touch.
        #------------------------------
        oidsmall = row[1]
        arcpy.SelectLayerByLocation_management(lyrSJ, 'SHARE_A_LINE_SEGMENT_WITH', row[0])
        break

#---------------------------------------------------------------
# Make sure the count is below the max threshold. If it is too
# high then remove touching features in descending count order.
#---------------------------------------------------------------
dictSelected = {}
with arcpy.da.SearchCursor(lyrSJ, ['OID@', 'Join_Count']) as cur:
    for row in cur:
        dictSelected[row[0]] = row[1]
count = sum(dictSelected.values())
if count > maxThreshold:
    while count > maxThreshold:
        s = sorted(dictSelected.items(), key=operator.itemgetter(1), reverse=True)
        del dictSelected[s[0][0]]
        count = sum(dictSelected.values())

    #-------------------------------
    # Reselect the features to use.
    #-------------------------------
    wcOID = 'OBJECTID in ('
    for k in dictSelected.keys():
        wcOID = wcOID + str(k) + ','
    wcOID = wcOID[0:-1] + ')'
    arcpy.MakeFeatureLayer_management(fcSJ, lyrSJ, wcOID)

#------------------------------
# Merge the features together.
#------------------------------
geomList = arcpy.Dissolve_management(lyrSJ, geom)

#-------------------------------
# Now update the feature class.
#-------------------------------
with arcpy.da.UpdateCursor(lyrSJ, flds) as cur:
    for row in cur:
        if row[1] == oidsmall:
            row[0] = geomList[0]
            row[2] = count
            cur.updateRow(row)
        else:
            cur.deleteRow()

Answer 2

I've outlined a process and given some sloppy code (I suggest you finish it, debug and post it as an answer when you are done). I think of this in three major steps, (1) iterate through the parcels (2) iterate through parcel's neighbors (3) split the last parcel if necessary. The third step is the hardest and I don't have a solution to that yet. Here is what I have so far:

import os
import arcpy
from arcpy import da
#da.SearchCursor is new as of 10.1 and is faster
parcels = "C:\Users\ln88615\Documents\GIS_data\cadastral_computation\mask2.shp"
parcels2 = path + "cpymask2.shp"
path = os.path.dirname(parcels)
parcelPoints = path + "tempPnts.shp"
outPnt = path + "tempPnt.shp"
pnts_layer = "points"    
parcel_layer = "parcel"
distanceTable = "distance"
#I've assumed this was your area field
area_field = "Matr"
thresholdArea = 1234
arcpy.FeatureToPoint_management(parcels_layer, parcelPoints)

The center points are to determine neighbors later. Select by location doesn't work as well since it will limit you to only those that border and you want to order your list as well, which distance is a good criteria for. Now you are ready to loop through the parcels

with arcpy.da.SearchCursor(parcels, ["FID", "area_field"]) as cursor:
    for row in cursor:
        area = row[1]
        if int(area) < thresholdArea:
            totalArea= int(area)
        #create layer for selected parcel
            sltqry = '"FID"=' + row[0]
            arcpy.MakeFeatureLayer_management(parcels, parcel_layer)
            arcpy.SelectLayerByAttribute_management(parcel_layer, "NEW_SELECTION", sltqry)
#(2) Within the above loop, you need to iterate through neighbor parcels.
            #Select point inside current parcel, must use layer
            arcpy.MakeFeatureLayer_management(parcelpoints, pnts_layer)
            arcpy.SelectLayerByLocation(pnts_layer, "WITHIN", parcel_layer)
        #It should do the select using the selection, if not create a FC after the select by attribute like:
        #arcpy.CopyFeatures_management(parcel_layer, path + "parcel.shp")
        #I created a new FC for the selection but PointDistance_analysis may take it more directly
            arcpy.CopyFeatures_management(pnts_layer, outPnt)
        #calculate distance to all other points, limit parcelPoints if it is long
            arcpy.PointDistance_analysis(outPnt, parcelPoints, distanceTable)
        #The first entry is its distance to itself so remove it or skip it in the next loop

You want to Join the distances back to the original parcels. Your table has IDs from the points which should still correspond to the ones from the parcels. Depending on how these are stored (it didn't look like you were using a gbd) you may or may not be able to rely on ID's matching. Look at your data and consider adding a "parcel_ID" to the points if need be.

            JoinField_management(parcels, "FID", distanceTable, "FID", "DISTANCE")
            #Sort the parcels by the distance, search cursor may do this
            arcpy.Sort_management(parcels, parcels2, "DISTANCE")
            neighborList = arcpy.da.SearchCursor(parcels2, ["DISTANCE", "area_field"])
            for neighbor in neighborList:
                if TotalArea<thresholdArea:
                    area2 = neighbor[1]
                    updateArea = totalArea + int(area2)
                    totalArea = updateArea
                    if totalArea>thresholdArea:
#                       (3) Do Split Function
            #outside the neighbor loop remove the join 
            arcpy.RemoveJoin_management (parcels)
        elseif int(area) > thresholdArea:
#           (3) Do Split Function