Explodierende Überlappung mit neuen nicht überlappenden Polygonen?

10

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?

ndimhypervol
quelle
Meinen Sie damit, dass Sie die Daten zu einem topologisch korrekten Satz "reduzieren" möchten?
Nagytech
@Geoist ZonalStats erfordert Polygone, die sich nicht überlappen. Wenn Sie eine überlappende Sammlung haben, besteht die offensichtliche, aber ineffiziente Lösung darin, die Polys zu durchlaufen und die Zonenstatistiken nacheinander zu berechnen. Es wäre effizienter, eine Teilmenge nicht überlappender Polys auszuwählen, Zonenstatistiken auf sie anzuwenden und zu iterieren. In der Frage wird gefragt, wie eine solche Auswahl effizient getroffen werden kann.
whuber
whuber - Ich denke, dass @Geoist vorschlägt, eine Reihe nicht überlappender Polygone aus den Schnittpunkten der Eingabepolygone zu erstellen. Schauen Sie sich dieses Bild an - (können Sie keine Bilder in Kommentaren posten?). Der Eingang befindet sich links. Die gesamte Region ist von drei Polygonen bedeckt, von denen jedes die beiden anderen schneidet. Die einzigen nicht überlappenden Teilmengen sind die Singletons, und diese erfüllen nicht Gotanukis Anforderung, dass die Gewerkschaft den Raum ausfüllt. Ich denke, Geoist schlägt vor, die Menge der sich nicht überschneidenden Regionen auf der rechten Seite zu erstellen, die für Zonenstatistiken gültig ist
Llaves
Ich denke, es gibt einige Verwirrung darüber, was das Endprodukt sein sollte. Könnten Sie ein Beispiel geben? Meine Interpretation ist, dass Sie möchten, dass die Ausgabe eine Auswahl von Polygonen ist, die sich nicht überlappen - während Sie die verbleibenden Polygone verwerfen oder auflösen. Arbeiten Sie mit einer oder mehreren Feature-Classes?
Aaron
1
Klingt für mich so, als ob @gotanuki die Mindestanzahl von Feature-Classes erstellen möchte, die nur nicht überlappende Polygone aus einer Polygon-Feature-Class mit überlappenden Polygonen enthalten
PolyGeo

Antworten:

14

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.

Karte von sieben 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

Vierfarbigkeit der Polygone

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.

whuber
quelle
Ich bin mir nicht sicher, ob dies die Frage beantwortet oder was die Frage ist, aber es ist trotzdem eine gute Antwort.
Nagytech
@Geoist Gibt es eine Möglichkeit, die Abbildung klarer zu machen oder das Problem besser zu erklären?
whuber
6

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
ndimhypervol
quelle
2
Um dies mit meiner Lösung zu verbinden: Ihre countOverlapsentspricht den beiden Zeilen nbrhoods <- sapply(vertices, neighbors); degrees <- sapply(nbrhoods, length)in meinem Code: degreesist die Überlappungszahl. Natürlich ist Ihr Code länger, da er den größten Teil der GIS-Analyse widerspiegelt, die in meiner Lösung als selbstverständlich angesehen wird: Sie identifizieren zunächst, welche Polygone sich überlappen, und am Ende verwenden Sie die Lösung zur Ausgabe von Polygon-Datasets. Es wäre eine gute Idee, die graphentheoretischen Berechnungen zu kapseln. Wenn Sie also jemals einen besseren
Farbalgorithmus
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)
Messamat
quelle
-3

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.

Alexandre Neto
quelle
2
Mit dieser Methode gelangen Sie nicht zum gewünschten Produkt. Hierbei handelt es sich um eine minimale Reihe von Auswahlen oder eindeutigen Funktionsklassen des Originals, die sich nicht überlappen. Die Produkte werden in die Zonenstatistik eingespeist. Daher ist es wichtig, die ursprüngliche Geometrie jedes Features beizubehalten.
Ndimhypervol
Du hast recht, sorry. Ich habe die Frage nicht gut verstanden. In diesem Fall und abhängig von der Größe des Rasters würde ich das Raster normalerweise in eine temporäre Punkt-Feature-Class (jede Zelle ein Punkt) konvertieren und eine räumliche Verknüpfung zwischen ihr und der Polygonebene durchführen. Vielleicht ist es ein sehr simpler und leistungsunfreundlicher Ansatz, aber er funktioniert und die überlappenden Polygone geben Ihnen kein Problem.
Alexandre Neto
Wenn ich richtig verstehe, was Sie mit dieser räumlichen Verknüpfung meinen, funktioniert Ihre zweite Lösung immer noch nicht, Alexandre, da zwischen den Punkten und den Polygonen eine Viele-zu-Viele-Beziehung besteht. Unabhängig davon ist dieser vektorbasierte Ansatz für jedes große Raster äußerst ineffizient und für große Raster unmöglich durchzuführen.
whuber
@whuber Sie haben Recht damit, ein sehr langsamer Prozess zu sein (Nehmen Sie mir eine halbe Stunde Zeit mit einem 4284 x 3009-Raster und 2401 Polygonen in einem Dualcore-RAM mit 2,8 GHz und 3 GB RAM mit Vista). Aber es funktioniert, da ich es schon getestet habe. In der räumlichen Verknüpfung müssen Sie eine Eins-zu-Eins-Beziehung verwenden und die Rasterwerte (als Mittelwert, Summe usw.) aggregieren. Das Ergebnis ist eine Vektorpolygon-Ebene ähnlich dem Original, jedoch mit einer neuen Spalte mit den aggregierten Rasterwerten, die jedes Polygon schneiden. Da dies keine optimale Lösung ist, kann dies für jemanden mit weniger Programmierkenntnissen nützlich sein (wie ich :-)).
Alexandre Neto