Explodierende Überlappung mit neuen nicht überlappenden Polygonen?

Angesichts mehrerer Polygone, die sich auf mehrere Arten überlappen, möchte ich aus diesen Features alle Polygone exportieren, die sich iterativ nicht mit anderen überlappen.

Das Produkt wäre eine Reihe von Merkmalen ohne Überlappung, die zusammengenommen das Original bilden.

Die Produkte könnten dann als Eingabe für die Zonenstatistik verwendet werden, und dies wäre viel schneller als das Iterieren der Zonenstatistik über jedes Polygon.

Ich habe versucht, dies in ArcPy ohne Erfolg zu codieren.

Gibt es dafür bereits Code?

Dies ist ein Problem beim Färben von Grafiken .

Denken Sie daran, dass eine Diagrammfärbung eine Zuordnung einer Farbe zu den Scheitelpunkten eines Diagramms ist, sodass keine zwei Scheitelpunkte, die eine Kante gemeinsam haben, dieselbe Farbe haben. Insbesondere sind die (abstrakten) Eckpunkte des Graphen die Polygone. Zwei Eckpunkte werden immer dann mit einer (ungerichteten) Kante verbunden, wenn sie sich schneiden (als Polygone). Wenn wir eine Lösung für das Problem finden - das eine Folge von (sagen wir k ) disjunkten Sammlungen der Polygone ist - und jeder Sammlung in der Folge eine eindeutige Farbe zuweisen, haben wir eine k- Färbung des Graphen erhalten . Es ist wünschenswert, ein kleines k zu finden .

Dieses Problem ist ziemlich schwierig und bleibt für beliebige Graphen ungelöst. Stellen Sie sich eine ungefähre Lösung vor, die einfach zu codieren ist. Ein sequentieller Algorithmus sollte reichen. Der Welsh-Powell-Algorithmus ist eine gierige Lösung, die auf einer absteigenden Reihenfolge der Eckpunkte nach Grad basiert. Sortiert in die Sprache der ursprünglichen Polygone, sortieren Sie die Polygone zunächst in absteigender Reihenfolge nach der Anzahl der anderen Polygone, die sie überlappen. Geben Sie dem ersten Polygon eine Anfangsfarbe. Versuchen Sie in jedem Schritt, das nächste Polygon mit einer vorhandenen Farbe zu färben. Wählen Sie also eine Farbe aus, die nicht vorhanden istbereits von einem der Nachbarn dieses Polygons verwendet. (Es gibt viele Möglichkeiten, zwischen den verfügbaren Farben zu wählen. Probieren Sie entweder die am wenigsten verwendete oder eine zufällig ausgewählte Farbe aus.) Wenn das nächste Polygon nicht mit einer vorhandenen Farbe gefärbt werden kann, erstellen Sie eine neue Farbe und färben Sie sie damit.

Wenn Sie eine Färbung mit einer kleinen Anzahl von Farben erreicht haben, führen Sie Zonenstatistiken Farbe für Farbe durch: Durch die Konstruktion wird garantiert, dass sich keine zwei Polygone einer bestimmten Farbe überlappen.

Hier ist Beispielcode in R. (Python-Code wäre nicht viel anders.) Zunächst beschreiben wir Überlappungen zwischen den sieben gezeigten Polygonen.

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)

Das heißt, die Polygone 1 und 2 überlappen sich, ebenso wie die Polygone 2 und 3, 3 und 4, ..., 1 und 7.

Sortieren Sie die Eckpunkte nach absteigendem Grad:

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

Ein (roher) sequentieller Farbalgorithmus verwendet die früheste verfügbare Farbe, die noch nicht von einem überlappenden Polygon verwendet wird:

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
}

Initialisieren Sie die Datenstrukturen ( colorsund color.next) und wenden Sie den Algorithmus an:

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

Teilen Sie die Polygone nach Farben in Gruppen auf:

split(vertices, colors)

Die Ausgabe in diesem Beispiel verwendet vier Farben:

$`1`
[1] 2 4

$`2`
[1] 3 6 7

$`3`
[1] 5

$`4`
[1] 1

Es hat die Polygone in vier nicht überlappende Gruppen unterteilt. In diesem Fall ist die Lösung nicht optimal ({{3,6,5}, {2,4}, {1,7}} ist für dieses Diagramm dreifarbig). Im Allgemeinen sollte die Lösung jedoch nicht schlecht sein.

Die von #whuber empfohlene Methodik hat mich dazu inspiriert, eine neue Richtung einzuschlagen, und hier ist meine bogenförmige Lösung in zwei Funktionen. Das erste, countOverlaps genannt, erstellt zwei Felder, "Overlaps" und "ovlpCount", um für jede Poly aufzuzeichnen, deren Polys sich überlappen und wie viele Überlappungen aufgetreten sind. Die zweite Funktion, explodeOverlaps, erstellt ein drittes Feld, "expl", das jeder Gruppe nicht überlappender Polys eine eindeutige Ganzzahl gibt. Der Benutzer kann dann basierend auf diesem Feld neue fc exportieren. Der Prozess ist in zwei Funktionen unterteilt, da sich das countOverlaps-Tool meiner Meinung nach als nützlich erweisen kann. Bitte entschuldigen Sie die Schlamperei des Codes (und die unachtsame Namenskonvention), da er ziemlich vorläufig ist, aber funktioniert. Stellen Sie außerdem sicher, dass der "idName" Feld stellt ein Feld mit eindeutigen IDs dar (nur mit ganzzahligen IDs getestet). Vielen Dank, dass Sie mir den Rahmen zur Verfügung gestellt haben, der zur Lösung dieses Problems erforderlich ist!

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

Es ist eine Weile her, aber ich habe diesen Code für meine eigene Anwendung verwendet und er hat großartig funktioniert - danke. Ich habe einen Teil davon neu geschrieben, um ihn zu aktualisieren, auf Zeilen anzuwenden (mit einer Toleranz) und ihn erheblich zu beschleunigen (siehe unten - ich verwende ihn auf 50 Millionen sich überschneidenden Puffern und es dauert nur ein paar Stunden).

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)

In diesen Fällen verwende ich im Allgemeinen die folgende Methode:

• Übergeben Sie die Feature-Class an eine UNION. (Es bricht die Polygone in all ihren Schnittpunkten)
• Fügen Sie die Felder X, Y und Area hinzu und berechnen Sie sie.
• Lösen Sie das Ergebnis nach X-, Y- und Flächenfeldern auf.

Ich glaube, das Ergebnis wird das sein, das Sie wollten, und Sie können sogar die Anzahl der Überlappungen zählen. Ich weiß nicht, ob es in Bezug auf die Leistung besser für Sie ist oder nicht.