1 votos

¿Fusión/disolución de todos los polígonos adyacentes entre sí mediante ArcGIS Desktop?

Tengo un shapefile que contiene muchos polígonos, que comparten límites o solo vértices, para polígonos diagonales.

Estoy tratando de fusionar/combinar todos los polígonos que comparten un límite.

He intentado usar Disolver seguido de un Multipart a Singlepart, sin embargo esto divide los polígonos que comparten un vértice (ver polígonos púrpura). Además, la parte superior izquierda de los polígonos púrpura tampoco se fusiona.

¿Hay alguna manera de hacer que, si un polígono toca el límite con otro polígono (utilizando el método D8, por lo que las 8 direcciones de una célula), debe ser fusionado.

enter image description here

1voto

Joel B Puntos 392

La extensión de Spatial Analyst ofrecería de nuevo una solución: utilizar el comando "Grupo de regiones" con el parámetro de conectividad de vecindad "OCHO" para asignar IDs únicos a cada zona rasterizada, y luego convertirla de nuevo a polígonos si es necesario.

Este flujo de trabajo continuaría en línea con su pregunta anterior Utilizar el buffer de IDW 3D Analyst para la propagación de la inundación

1voto

UnkwnTech Puntos 21942

No voy a decir que entiendo cada parte de mi solución porque una gran parte del código en el medio viene de un Q&A que inicié en Stack Overflow: https://stackoverflow.com/questions/49396943

El código de la solución comienza creando algunos datos de prueba que se ilustran etiquetando sus polígonos con sus OBJECTIDs:

enter image description here

import arcpy

# Create test data
if arcpy.Exists("C:/Temp/test.gdb"):
    arcpy.Delete_management("C:/Temp/test.gdb")
arcpy.CreateFileGDB_management("C:/Temp","test")
arcpy.CreateFishnet_management("C:/Temp/test.gdb/testFishnet",
                               "0 0","0 1","1","1","5","5",
                               labels="NO_LABELS",
                               geometry_type="POLYGON")
arcpy.Select_analysis("C:/Temp/test.gdb/testFishnet",
                      "C:/Temp/test.gdb/testFC",
                      "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
arcpy.AddField_management("C:/Temp/test.gdb/testFC","DISS","LONG")

# Write neighbourhood relationships to a table
arcpy.PolygonNeighbors_analysis("C:/Temp/test.gdb/testFC",
                                "C:/Temp/test.gdb/testNeighbors",
                                "OBJECTID","NO_AREA_OVERLAP",
                                "BOTH_SIDES")

# 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]]
        else:
            d[row[0]] = [row[1]]
print(d)
# 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
    else:
        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
    else:
        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]:
        edges.append((u,vertices[v]))

# do disjoint set union find similar to kruskal's algorithm
for u,v in edges:
    if find_set(u) != find_set(v):
        union(u,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)

print(parents)
# 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]]
        cursor.updateRow(row)
# Dissolve on the dissolve field (DISS)
arcpy.Dissolve_management("C:/Temp/test.gdb/testFC",
                          "C:/Temp/test.gdb/testDiss",
                          "DISS","","MULTI_PART","DISSOLVE_LINES")

El resultado de la ejecución del código es una clase de característica disuelta en los identificadores de grupo que asigna en función de si las características de entrada se tocan por al menos un nodo.

enter image description here

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X