¿Explosión superpuesta a nuevos polígonos no superpuestos?

Dados múltiples polígonos que se superponen de múltiples maneras, me gustaría exportar de estas características todos los polígonos que no se superponen con otros, de forma iterativa.

El producto sería una serie de características sin solapamiento que, cuando se suman, forman el original.

Los productos podrían utilizarse como entrada para las estadísticas zonales, y esto sería mucho más rápido que iterar estadísticas zonales sobre cada polígono.

He estado tratando de codificar esto en ArcPy sin éxito.

¿El código para hacer esto ya existe?

Este es un problema de coloración gráfica .

Recuerde que un color de gráfico es una asignación de un color a los vértices de un gráfico de tal manera que no haya dos vértices que compartan un borde también tendrán el mismo color. Específicamente, los vértices (abstractos) del gráfico son los polígonos. Dos vértices están conectados con un borde (no dirigido) cada vez que se cruzan (como polígonos). Si tomamos alguna solución al problema, que es una secuencia de (digamos k ) colecciones disjuntas de los polígonos, y asignamos un color único a cada colección en la secuencia, entonces habremos obtenido una k- coloración del gráfico . Es deseable encontrar una pequeña k .

Este problema es bastante difícil y sigue sin resolverse para gráficos arbitrarios. Considere una solución aproximada que sea simple de codificar. Un algoritmo secuencial debe hacer. El algoritmo galés-Powell es una solución codiciosa basada en un orden descendente de los vértices por grado. Traducido al lenguaje de los polígonos originales, primero ordene los polígonos en orden descendente del número de otros polígonos que se superponen. Trabajando en orden, déle al primer polígono un color inicial. En cada paso sucesivo, intente colorear el siguiente polígono con un color existente: es decir, elija un color que no seaya utilizado por cualquiera de los vecinos de ese polígono. (Hay muchas formas de elegir entre los colores disponibles; pruebe el que menos se haya utilizado o elija uno al azar). Si el próximo polígono no se puede colorear con un color existente, cree un nuevo color y asígnelo.

Una vez que haya logrado una coloración con una pequeña cantidad de colores, realice las zonas zonal color por color: por construcción, se garantiza que no se superponen dos polígonos de un color determinado.

Aquí hay un código de muestra R. (El código de Python no sería muy diferente). Primero, describimos las superposiciones entre los siete polígonos mostrados.

edges <- matrix(c(1,2, 2,3, 3,4, 4,5, 5,1, 2,6, 4,6, 4,7, 5,7, 1,7), ncol=2, byrow=TRUE)

Es decir, los polígonos 1 y 2 se superponen, y también lo hacen los polígonos 2 y 3, 3 y 4, ..., 1 y 7.

Ordenar los vértices por grado descendente:

vertices <- unique(as.vector(edges))
neighbors <- function(i) union(edges[edges[, 1]==i,2], edges[edges[, 2]==i,1])
nbrhoods <- sapply(vertices, neighbors)
degrees <- sapply(nbrhoods, length)
v <- vertices[rev(order(degrees))]

Un algoritmo de coloración secuencial (en bruto) utiliza el primer color disponible que no haya sido utilizado por ningún polígono superpuesto:

color <- function(i) {
  n <- neighbors(i)
  candidate <- min(setdiff(1:color.next, colors[n]))
  if (candidate==color.next) color.next <<- color.next+1
  colors[i] <<- candidate
}

Inicialice las estructuras de datos ( colorsy color.next) y aplique el algoritmo:

colors <- rep(0, length(vertices))
color.next <- 1
temp <- sapply(v, color)

Divide los polígonos en grupos según el color:

split(vertices, colors)

La salida en este ejemplo usa cuatro colores:

$`1`
[1] 2 4

$`2`
[1] 3 6 7

$`3`
[1] 5

$`4`
[1] 1

Ha dividido los polígonos en cuatro grupos no superpuestos. En este caso, la solución no es óptima ({{3,6,5}, {2,4}, {1,7}} son tres colores para este gráfico). Sin embargo, en general, la solución que obtiene no debería ser tan mala.

La metodología recomendada por #whuber me inspiró a tomar una nueva dirección, y aquí está mi solución arcpy, en dos funciones. El primero, llamado countOverlaps, crea dos campos, "overlaps" y "ovlpCount" para registrar para cada poli que polys se superpuso con él y cuántas superposiciones ocurrieron. La segunda función, explodeOverlaps, crea un tercer campo, "expl", que proporciona un número entero único a cada grupo de polys no superpuestos. El usuario puede exportar nuevos fc basados ​​en este campo. El proceso se divide en dos funciones porque creo que la herramienta countOverlaps puede resultar útil por sí misma. Disculpe la descuido del código (y la convención de nomenclatura descuidada), ya que es bastante preliminar, pero funciona. También asegúrese de que el "idName" El campo representa un campo con ID únicos (solo probado con ID enteros). ¡Gracias whuber por proporcionarme el marco necesario para abordar este problema!

def countOverlaps(fc,idName):
    intersect = arcpy.Intersect_analysis(fc,'intersect')
    findID = arcpy.FindIdentical_management(intersect,"explFindID","Shape")
    arcpy.MakeFeatureLayer_management(intersect,"intlyr")
    arcpy.AddJoin_management("intlyr",arcpy.Describe("intlyr").OIDfieldName,findID,"IN_FID","KEEP_ALL")
    segIDs = {}
    featseqName = "explFindID.FEAT_SEQ"
    idNewName = "intersect."+idName

    for row in arcpy.SearchCursor("intlyr"):
        idVal = row.getValue(idNewName)
        featseqVal = row.getValue(featseqName)
        segIDs[featseqVal] = []
    for row in arcpy.SearchCursor("intlyr"):
        idVal = row.getValue(idNewName)
        featseqVal = row.getValue(featseqName)
        segIDs[featseqVal].append(idVal)

    segIDs2 = {}
    for row in arcpy.SearchCursor("intlyr"):
        idVal = row.getValue(idNewName)
        segIDs2[idVal] = []

    for x,y in segIDs.iteritems():
        for segID in y:
            segIDs2[segID].extend([k for k in y if k != segID])

    for x,y in segIDs2.iteritems():
        segIDs2[x] = list(set(y))

    arcpy.RemoveJoin_management("intlyr",arcpy.Describe(findID).name)

    if 'overlaps' not in [k.name for k in arcpy.ListFields(fc)]:
        arcpy.AddField_management(fc,'overlaps',"TEXT")
    if 'ovlpCount' not in [k.name for k in arcpy.ListFields(fc)]:
        arcpy.AddField_management(fc,'ovlpCount',"SHORT")

    urows = arcpy.UpdateCursor(fc)
    for urow in urows:
        idVal = urow.getValue(idName)
        if segIDs2.get(idVal):
            urow.overlaps = str(segIDs2[idVal]).strip('[]')
            urow.ovlpCount = len(segIDs2[idVal])
        urows.updateRow(urow)

def explodeOverlaps(fc,idName):

    countOverlaps(fc,idName)

    arcpy.AddField_management(fc,'expl',"SHORT")

    urows = arcpy.UpdateCursor(fc,'"overlaps" IS NULL')
    for urow in urows:
        urow.expl = 1
        urows.updateRow(urow)

    i=1
    lyr = arcpy.MakeFeatureLayer_management(fc)
    while int(arcpy.GetCount_management(arcpy.SelectLayerByAttribute_management(lyr,"NEW_SELECTION",'"expl" IS NULL')).getOutput(0)) > 0:
        ovList=[]
        urows = arcpy.UpdateCursor(fc,'"expl" IS NULL','','','ovlpCount D')
        for urow in urows:
            ovVal = urow.overlaps
            idVal = urow.getValue(idName)
            intList = ovVal.replace(' ','').split(',')
            for x in intList:
                intList[intList.index(x)] = int(x)
            if idVal not in ovList:
                urow.expl = i
            urows.updateRow(urow)
            ovList.extend(intList)
        i+=1

Ha pasado un tiempo, pero usé este código para mi propia aplicación y ha funcionado muy bien, gracias. Reescribí parte de él para actualizarlo, aplicarlo a las líneas (con tolerancia) y acelerarlo significativamente (a continuación, lo estoy ejecutando en 50 millones de búferes de intersección y solo toma un par de horas).

def ExplodeOverlappingLines(fc, tolerance, keep=True):
        print('Buffering lines...')
        idName = "ORIG_FID"
        fcbuf = arcpy.Buffer_analysis(fc, fc+'buf', tolerance, line_side='FULL', line_end_type='FLAT')
        print('Intersecting buffers...')
        intersect = arcpy.Intersect_analysis(fcbuf,'intersect')

        print('Creating dictionary of overlaps...')
        #Find identical shapes and put them together in a dictionary, unique shapes will only have one value
        segIDs = defaultdict(list)
        with arcpy.da.SearchCursor(intersect, ['Shape@WKT', idName]) as cursor:
            x=0
            for row in cursor:
                if x%100000 == 0:
                    print('Processed {} records for duplicate shapes...'.format(x))
                segIDs[row[0]].append(row[1])
                x+=1

        #Build dictionary of all buffers overlapping each buffer
        segIDs2 = defaultdict(list)
        for v in segIDs.values():
            for segID in v:
                segIDs2[segID].extend([k for k in v if k != segID and k not in segIDs2[segID]])

        print('Assigning lines to non-overlapping sets...')
        grpdict = {}
        # Mark all non-overlapping one to group 1
        for row in arcpy.da.SearchCursor(fcbuf, [idName]):
            if row[0] in segIDs2:
                grpdict[row[0]] = None
            else:
                grpdict[row[0]] = 1

        segIDs2sort = sorted(segIDs2.items(), key=lambda x: (len(x[1]), x[0])) #Sort dictionary by number of overlapping features then by keys
        i = 2
        while None in grpdict.values(): #As long as there remain features not assigned to a group
            print(i)
            ovset = set()  # list of all features overlapping features within current group
            s_update = ovset.update
            for rec in segIDs2sort:
                if grpdict[rec[0]] is None: #If feature has not been assigned a group
                    if rec[0] not in ovset: #If does not overlap with a feature in that group
                        grpdict[rec[0]] = i  # Assign current group to feature
                        s_update(rec[1])  # Add all overlapping feature to ovList
            i += 1 #Iterate to the next group

        print('Writing out results to "expl" field in...'.format(fc))
        arcpy.AddField_management(fc, 'expl', "SHORT")
        with arcpy.da.UpdateCursor(fc,
                                   [arcpy.Describe(fc).OIDfieldName, 'expl']) as cursor:
            for row in cursor:
                if row[0] in grpdict:
                    row[1] = grpdict[row[0]]
                    cursor.updateRow(row)

        if keep == False:
            print('Deleting intermediate outputs...')
            for fc in ['intersect', "explFindID"]:
                arcpy.Delete_management(fc)

En estos casos, generalmente uso el siguiente método:

• Pase la clase de entidad a través de una UNIÓN; (Rompe los polígonos en todas sus intersecciones)
• Agregue los campos X, Y y Área y calcule;
• Disuelva el resultado por X, Y, campos de área.

Creo que el resultado será el que deseabas, e incluso puedes contar el número de superposiciones. No sé si en términos de rendimiento será mejor para usted o no.