Eksploduje nakładanie się na nowe nie nakładające się wielokąty?

Biorąc pod uwagę wiele wielokątów, które nakładają się na wiele sposobów, chciałbym wyeksportować z tych funkcji wszystkie wielokąty, które nie pokrywają się iteracyjnie.

Produkt miałby szereg funkcji bez nakładania się, które po zsumowaniu składają się na oryginał.

Produkty mogłyby następnie zostać wykorzystane jako dane wejściowe do statystyki strefowej, a byłoby to znacznie szybsze niż iteracja statystyki strefowej dla każdego wielokąta.

Próbowałem to zakodować w ArcPy bez powodzenia.

Czy kod do tego celu już istnieje?

Jest to problem z kolorowaniem wykresów .

Przypomnij sobie, że kolorowanie wykresu to przypisanie koloru do wierzchołków wykresu w taki sposób, że żadne dwa wierzchołki, które dzielą krawędź, również nie będą miały tego samego koloru. W szczególności (abstrakcyjne) wierzchołki wykresu to wielokąty. Dwa wierzchołki są połączone (niekierowaną) krawędzią, gdy przecinają się (jako wielokąty). Jeśli weźmiemy jakieś rozwiązanie problemu - który jest sekwencją (powiedzmy k ) rozłącznych kolekcji wielokątów - i przypiszemy unikalny kolor do każdej kolekcji w sekwencji, uzyskamy k -kolorowanie wykresu . Pożądane jest znalezienie małego k .

Ten problem jest dość trudny i pozostaje nierozwiązany w przypadku dowolnych grafów. Rozważ przybliżone rozwiązanie, które można łatwo kodować. Algorytm sekwencyjny powinien zrobić. Algorytm Welsha-Powella jest chciwym rozwiązaniem opartym na malejącej kolejności wierzchołków według stopnia. Przetłumaczone na język oryginalnych wielokątów, najpierw posortuj wielokąty w porządku malejącym względem liczby innych wielokątów, które się na siebie nakładają. Pracując w porządku, nadaj pierwszemu wielokątowi początkowy kolor. Na każdym kolejnym kroku spróbuj pokolorować następny wielokąt istniejącym kolorem: wybierz kolor, który nie jestużywany już przez któregokolwiek z sąsiadów tego wielokąta. (Istnieje wiele sposobów wyboru spośród dostępnych kolorów; wypróbuj ten, który był najmniej używany, lub wybierz jeden losowo.) Jeśli następnego wielokąta nie można pokolorować istniejącym kolorem, utwórz nowy kolor i pokoloruj go tym.

Po uzyskaniu kolorowania za pomocą niewielkiej liczby kolorów wykonaj statystyki strefowe kolor po kolorze: z uwagi na konstrukcję masz pewność, że żadne dwa wielokąty danego koloru nie zachodzą na siebie.

Oto przykładowy kod w R. (Kod w Pythonie nie różni się zbytnio.) Po pierwsze, opisujemy nakładanie się siedmiu pokazanych wielokątów.

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)

Oznacza to, że wielokąty 1 i 2 nakładają się, podobnie jak wielokąty 2 i 3, 3 i 4, ..., 1 i 7.

Sortuj wierzchołki według malejącego stopnia:

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

(Surowy) algorytm sekwencyjnego kolorowania używa najwcześniejszego dostępnego koloru, który nie jest jeszcze używany przez nakładający się wielokąt:

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
}

Zainicjuj struktury danych ( colorsi color.next) i zastosuj algorytm:

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

Podziel wielokąty na grupy według kolorów:

split(vertices, colors)

Dane wyjściowe w tym przykładzie wykorzystują cztery kolory:

$`1`
[1] 2 4

$`2`
[1] 3 6 7

$`3`
[1] 5

$`4`
[1] 1

Podzielił wielokąty na cztery niezachodzące na siebie grupy. W tym przypadku rozwiązanie nie jest optymalne ({{3,6,5}, {2,4}, {1,7}} to trójkolorowy wykres). Ogólnie rzecz biorąc, rozwiązanie, które dostaje, nie powinno być jednak takie złe.

Metodologia zalecana przez #whuber zainspirowała mnie do obrania nowego kierunku, a oto moje arkadowe rozwiązanie w dwóch funkcjach. Pierwsze, zwane countOverlaps, tworzy dwa pola: „overlaps” i „ovlpCount”, aby zapisać dla każdego poli, które pokrywały się z nim, i ile wystąpiło. Druga funkcja, explodeOverlaps, tworzy trzecie pole, „expl”, które daje unikalną liczbę całkowitą dla każdej grupy nie nakładających się polis. Użytkownik może następnie wyeksportować nowe fc na podstawie tego pola. Proces jest podzielony na dwie funkcje, ponieważ myślę, że narzędzie countOverlaps może się przydać. Proszę wybaczyć niechlujność kodu (i nieostrożną konwencję nazewnictwa), ponieważ jest to dość wstępne, ale działa. Upewnij się także, że „idName” pole reprezentuje pole z unikalnymi identyfikatorami (testowane tylko z identyfikatorami liczb całkowitych). Dziękuję fanom za udostępnienie mi ram niezbędnych do rozwiązania tego problemu!

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

Minęło trochę czasu, ale użyłem tego kodu do własnej aplikacji i działa świetnie - dziękuję. Ponownie napisałem część, aby ją zaktualizować, zastosować do linii (z tolerancją) i znacznie przyspieszyć (poniżej - korzystam z 50 milionów przecinających się buforów i zajmuje to tylko kilka godzin).

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)

W takich przypadkach zazwyczaj używam następującej metody:

• Przekaż klasę obiektów przez UNION; (Łamie wielokąty na wszystkich skrzyżowaniach)
• Dodaj pola X, Y i Obszar i oblicz je;
• Rozpuść wynik za pomocą pól X, Y, Obszar.

Wierzę, że wynik będzie taki, jaki chciałeś, i możesz nawet policzyć liczbę nakładek. Nie wiem, czy pod względem wydajności będzie to dla Ciebie lepsze, czy nie.