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.enter image description here

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")
  • It seems like there will be several points in this process you will need to decide what you want the script to do. First off, what order should it go through the neighboring parcels, 520 or 692 first? Then when you split the last parcel where should it put the split?
    – jbosq
    Commented Nov 20, 2014 at 16:05
  • Thanks for your comments. The order of searching neighbor could be any direction , just making sure that all the neighbor parcels (having common boundary) gets selected first. For example, 520 and 692 should be selected first instead of the order 692 and 450. Splitting polygons can be done by selecting the middle point of the parcel. then shifting it in an iterative manner. Thanks
    – Zaman
    Commented Nov 20, 2014 at 17:04
  • What version of arc are you using (10.0 / 10.2)? I'm assuming you don't have code already and are starting from scratch?
    – jbosq
    Commented Nov 20, 2014 at 17:49
  • I have arcGIS 10.2 and want to start from scratch. thanks
    – Zaman
    Commented Nov 20, 2014 at 17:51
  • Now I have added some preliminary code. This python script can select all the neighboring parcels for a given selected one. But it is still not region-growing. For example, after selecting parcel 22 (neighbors are 692, 520) , when I select parcel 592 (neighbors are 22, 692 and overlaps with 692 from previous one). How can I still select one by one parcel. Thanks
    – Zaman
    Commented Nov 21, 2014 at 9:41

2 Answers 2


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])

# 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

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
  • I got stuck in the second part where point distance is implemented. What is the difference between outPnt and parcelPoints? What is the function of sltqry?
    – Zaman
    Commented Nov 26, 2014 at 13:34
  • parcelPoints has a point for each parcel. You want outPnt to be just the point corresponding to the parcel currently being processed. This way point distance does the distance to your point of interest as it moves through the parcels rather than having to calculate distances from all points to all points. sltqry is the where clause for a selection based on the FID. Try the edits I made above. Alternatively you may be able to use sltqry to select the point you want and create outPnt more directly as long as point 1 of parcelPoints is located inside parcel 1.
    – jbosq
    Commented Nov 26, 2014 at 20:40
  • The PointDistance_analysis is not workign for some unknown reason and the following error shows up:ERROR 000210: Cannot create output distance Failed to execute (PointDistance).
    – Zaman
    Commented Dec 12, 2014 at 15:18
  • Whoops, the file path needed an extension (distanceTable = "distance.txt"). If that doesn't fix it outPnt and parcelPoints should be created as FCs before it fails, so check that PointDistance works manually outside the code with your same inputs. Also, this question should help you do the split function
    – jbosq
    Commented Dec 19, 2014 at 17:54
  • I got neighborList implemented and got totalArea as desired. But how do I select parcels sequentially when the algorithm stops just before threshold value? Do I have to implement SelectLayer with attribute?
    – Zaman
    Commented Jan 7, 2015 at 16:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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