I have a shapefile containing a lot of polygons, that share boundaries or just single vertices, for diagonal polygons.

I'm trying to merge/combine all the polygons that share a boundary.

I've tried using Dissolve followed by a Multipart to Singlepart, however this splits up polygons that share a vertice (see purple polygons). Also the top left part of the purple polygons aren't merged either.

Is there a way to make it so, if a polygon touches boundary with another polygon (using D8 method, so all 8 directions from a cell), it should be merged.

enter image description here

  • Couldn't you make a small buffer of your polygon features that you dissolve then multi to single then spatial join you buffer fid to your polygons and merge the polygons based on that id? – Boodoo Mar 16 '18 at 9:12
  • You could try Integrate and then dissolve – BERA Mar 16 '18 at 9:13
  • 1
    Quick thought is to get your neighborhood relationships into a table using Polygon Neighbors and then use search and update cursors to read it and set a field to values you can dissolve on in the way you want. Polygon Neighbors knows about both node and arc neighbors. – PolyGeo Mar 16 '18 at 9:21
  • @PolyGeo This sounds about right. Though I'm not quite sure how to do the Search and update you describe, though I feel like it should be fairly logical I can't quite grasp it. When I've got the src_ObjectID and nbr_objectID, how do i proceed with the updating? – FoolzRailer Mar 20 '18 at 8:19
  • I'm still looking into this but have hit a gap in my Python skills that I am try to get past at stackoverflow.com/questions/49396943/… – PolyGeo Mar 21 '18 at 2:09

Spatial Analyst extension would again offer a solution: use the “Region Group” command with the “EIGHT” neighbourhood connectivity parameter to assign unique IDs this each raster zone, then convert back to polygons if needed.

This workflow would continue in line with your previous question Use buffer from IDW 3D Analyst for flood spread


I won't claim to understand every part of my solution because a large chunk of code in the middle comes from a Q&A that I initiated at Stack Overflow: https://stackoverflow.com/questions/49396943

The solution code starts by creating some test data which is illustrated by labelling its polygons with their OBJECTIDs:

enter image description here

import arcpy

# Create test data
if arcpy.Exists("C:/Temp/test.gdb"):
                               "0 0","0 1","1","1","5","5",
                      "OID IN (1,7,13,19,25,16,21,22,4,5,10)")

# Add field that will hold group numbers and will eventually be dissolved on

# Write neighbourhood relationships to a table

# Write a dictionary with src_OBJECTID as its keys and a list of
# nbr_OBJECTIDs as its values
d = {}
with arcpy.da.SearchCursor("C:/Temp/test.gdb/testNeighbors",
                           ["src_OBJECTID","nbr_OBJECTID"]) as cursor:
    for row in cursor:
        if row[0] in d:
            d[row[0]] = d[row[0]] + [row[1]]
            d[row[0]] = [row[1]]
# For the test data the dictionary (d) looks like this
# {1: [4], 2: [3, 5], 3: [2, 5], 4: [1, 6], 5: [2, 3], 
#  6: [4, 8], 7: [9, 10], 8: [6, 11], 9: [7, 10], 10: [7, 9], 11: [8]}

# Code below here comes from https://stackoverflow.com/a/49397534/820534
# with minor modification because naming groups with integers rather
# letters is easier
class MSet(object):
    def __init__(self, p):
        self.val = p
        self.p = self
        self.rank = 0

def parent_of(x): # recursively find the parents of x
    if x.p == x:
        return x.val
        return parent_of(x.p)

def make_set(x):
    return MSet(x)

def find_set(x):
    if x != x.p:
        x.p = find_set(x.p)
    return x.p

def link(x,y):
    if x.rank > y.rank:
        y.p = x
        x.p = y
        if x.rank == y.rank:
            y.rank += 1

def union(x,y):
    link(find_set(x), find_set(y))

vertices = {k: make_set(k) for k in d.keys()}
edges = []

for k,u in vertices.items():
    for v in d[k]:

# do disjoint set union find similar to kruskal's algorithm
for u,v in edges:
    if find_set(u) != find_set(v):

# resolve the root of each disjoint set
parents = {}

# generate set of parents
set_parents = set() 

for u,v in edges:
    set_parents |= {parent_of(u)}
    set_parents |= {parent_of(v)}

# make a mapping from only parents to A-Z, to allow up to 26 disjoint sets
# letters = {k : chr(v) for k,v in zip(set_parents, list(range(65,91)))}

for u,v in edges:
#     parents[u.val] = letters[parent_of(u)]
#     parents[v.val] = letters[parent_of(v)]
    parents[u.val] = parent_of(u)
    parents[v.val] = parent_of(v)

# End of code from https://stackoverflow.com/a/49397534/820534

# For the test data the dictionary (parents) looks like this
# {1: 4, 2: 3, 3: 3, 4: 4, 5: 3, 6: 4, 7: 9, 8: 4, 9: 9, 10: 9, 11: 4}

# Read the parents dictionary and update the dissolve field
# in the original feature class
with arcpy.da.UpdateCursor("C:/Temp/test.gdb/testFC",
                           ["OBJECTID","DISS"]) as cursor:
    for row in cursor:
        row[1] = parents[row[0]]
# Dissolve on the dissolve field (DISS)

The result of running the code is a feature class dissolved on the group identifiers that it assigns based on whether input features touch by at least a node.

enter image description here

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.