4

I need to numerate a thousands of polygons in groups, so each group will get unique number. Every intersecting polygons is unique group (see picture).

So polygons have field "groups" with values:

polygons 1,2,3 is group 1;

polygons 4,5 is group 2;

polygons 6,7,8,9 is group 3 and so on...

enter image description here

My current approach is very slow it takes about 5 minutes to enumerate layer with about 10 000 polygons (it is 4 100 groups).

Is there any faster way to do this?

Here is what I do:

  1. Dissolve all polygons - that is my way to find groups (arcpy.Dissolve_management)
  2. Give unique numbers to every polygon in dissolved layer (arcpy.CalculateField_management)
  3. Iterate every feature in dissolved layer (arcpy.SelectLayerByAttribute_management)
  4. Select polygons that intersect with original (not dissolved) layer.
  5. Give them group number (arcpy.CalculateField_management)

Here is a part of the code:

#1
arcpy.Dissolve_management(fc, dissolve_shp, "", "", "SINGLE_PART")

#2
arcpy.CalculateField_management(dissolve_shp, "ID", '!FID!+1', "PYTHON")

#3
number_of_rows = arcpy.GetCount_management(dissolve_shp)
number_of_rows = int(number_of_rows.getOutput(0))
lyr_dissolve = r"in_memory\Temporary_layer_DISSOLVE_DEL"
lyr_merge = r"in_memory\Temporary_layer_MERGE_DEL"    
arcpy.MakeFeatureLayer_management(dissolve_shp, lyr_dissolve)
arcpy.MakeFeatureLayer_management(fc, lyr_merge)

for i in range(number_of_rows):
    i2=i+1
    if float(i2)//100 == float(i2)/100:
        print i2, "of", number_of_rows

    arcpy.SelectLayerByAttribute_management (lyr_dissolve, "NEW_SELECTION", '"ID" = '+str(i2))

    #4
    arcpy.SelectLayerByLocation_management (lyr_merge, "INTERSECT", lyr_dissolve)

    #5
    arcpy.CalculateField_management(lyr_merge, "TERM_P_NUM", i2, "PYTHON")
  • 5
    Seeing as the groups are the singles dissolved why not Intersect (or union) them, that gets the GroupID onto the polygons... adding a field (OrigID) and copying FID values prior to union will mean you can find the original (source) polygon for the intersection area. Then iterate the polygons with an UpdateCursor using whereclause GroupID = 1..n to enumerate each polygon in the group. – Michael Stimson Jun 17 '15 at 5:31
  • 1
    @MichaelMiles-Stimson I haven't run a test to confirm but I think that polygon overlay should allow the asker to dispense with the ArcPy cursor and thousands of Select Layer By Locations, to bring the time down to around 10-15 seconds in total. They owe you one. I don't think they are asking to do that last enumeration that you mention. – PolyGeo Jun 17 '15 at 5:37
  • Can you access to Polygon Neighbors tool under Analysis>Proximity toolbox? – fatih_dur Jun 17 '15 at 6:02
  • @ fatih_dur yes i have all ArcGIS tools available – Mr. Che Jun 17 '15 at 6:08
  • @MichaelMiles-Stimson Your approach works fine and it is much faster than mine. The only problem is splitting polygons. Union tool splits all polygons by intersection line. I bypass this by dissolving it back (use attribute fields option). – Mr. Che Jun 17 '15 at 12:10
2

A faster way to accomplish this is with a spatial join for step 3, once your polygons have their unique number. Let's say UID is the field name for your unique ID from your dissolve polygons, Dissolve is your dissolve polygons, Polies is your target polygon, and outFC will be the full path of your resulting feature class.

Below I get a little fancy with field maps to limit the fields returned. Note I haven't tested the script.

#Create field mappings object
fms = arcpy.FieldMappings()
#Add fields from target
fms.addTable (Polies)

#Add UID field from dissolve
uidFm = arcpy.FieldMap ()
uidFm.addInputField (Dissolve, UID)
fms.addFieldMap (uidFm)

#Spatial join
arcpy.SpatialJoin_analysis (Polies, Dissolve, outFC, field_mapping = fms) 
  • +1 A good approach. I would love to see a benchmark to see how this method compares to the OP's;) – Aaron Jun 18 '15 at 2:53
1

Here is another approach by using PolygonNeighbors_analysis, which was introduced at ArcGIS 10.1.

# fc                : Feature class to-be-processed
# ID_field          : Unique ID field in the feature class on which groupping will be based
# neighbour_table   : Table location where the Polygon Neighbor results will be saved,
#                     preferrably a GDB location (considering field naming limitation of dbf)
arcpy.PolygonNeighbors_analysis(fc,neighbour_table,ID_field,"NO_AREA_OVERLAP","BOTH_SIDES","#","METERS","SQUARE_METERS")
group_ID=0
pairs=[]
for row in arcpy.SearchCursor(neighbour_table):
    pairs.append([row.getValue('src_'+ID_field),row.getValue('nbr_'+ID_field)])

pairs=sorted(pairs)
group_dict={group_ID:set(pairs[0])}
for pair in pairs[1:]:
    scr,nbr=pair
    processed=False
    for k in group_dict.keys():
        if not processed:
            if scr in group_dict[k] or nbr in group_dict[k]:
                group_dict[k]=group_dict[k].union([scr,nbr])
                processed=True
        else:
            break
    if not processed:
        group_ID+=1
        group_dict[group_ID]=set([scr,nbr])

group_dict_invert=dict((val,k) for k in group_dict.keys() for val in group_dict[k])
groupping_field="Group_ID"
arcpy.AddField_management(fc,groupping_field,"Long")
up_cur=arcpy.UpdateCursor(fc)
for row in up_cur:
    ID=row.getValue(ID_field)
    if ID in group_dict_invert.keys():
        row.setValue(groupping_field,group_dict_invert[ID])
    else:
        group_ID+=1
        row.setValue(groupping_field,group_ID)
    up_cur.updateRow(row)

del row,up_cur

enter image description here

  • +1 You will see speed improvements if you use the da cursors. It would be interesting to see a benchmark test using this method vs the OP's;) – Aaron Jun 18 '15 at 2:50
  • 1
    Polygon Neighbors does not need Advanced - for at least last few versions Basic is sufficient and I think since tool was introduced at about 10.1. – PolyGeo Jun 18 '15 at 11:10
  • 1
    I quickly edited a screenshot (from the link you provided) into your answer to illustrate that it is licensed at all levels. – PolyGeo Jun 18 '15 at 23:16

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.