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

10

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?

ndimhypervol
źródło
Czy masz na myśli, że chcesz „spłaszczyć” dane do zestawu poprawnego topologicznie?
nagytech
@Geoist ZonalStats wymaga wielokątów, które się nie nakładają. Kiedy masz nakładającą się kolekcję, oczywistym, ale nieefektywnym rozwiązaniem jest zapętlenie polis i obliczenie statystyk strefowych jeden po drugim. Bardziej efektywne byłoby wybranie podzbioru nienakładających się polis, zastosowanie do nich statystyk strefowych i iteracja. Pytanie dotyczy tego, jak skutecznie dokonywać takich wyborów.
whuber
whuber - Myślę, że @Geoist sugeruje utworzenie zestawu nienakładających się wielokątów ze skrzyżowań wielokątów wejściowych. Spójrz na to zdjęcie - (nie możesz dodawać zdjęć w komentarzach?). Wejście znajduje się po lewej stronie. Cały region jest pokryty trzema wielokątami, z których każdy przecina oba pozostałe. Jedynymi nienakładającymi się podzbiorami są singletony, które nie spełniają wymogu Gotanuki, że związek wypełnia przestrzeń. Myślę, że Geoist sugeruje utworzenie zestawu nie przecinających się regionów po prawej stronie, który byłby ważny dla stref strefowych
Llaves
Myślę, że istnieje pewne zamieszanie co do tego, jaki powinien być produkt końcowy. Czy możesz podać przykład? Moja interpretacja jest taka, że ​​chciałbyś, aby wynik był zbiorem wielokątów, które nie zachodzą na siebie - podczas odrzucania lub rozpuszczania pozostałych wielokątów. Czy pracujesz z jedną lub wieloma klasami funkcji?
Aaron
1
Wydaje mi się, że @gotanuki chce stworzyć minimalną liczbę klas obiektów, które zawierają tylko nie nakładające się wielokąty z klasy obiektów wieloboków z nakładającymi się wielokątami
PolyGeo

Odpowiedzi:

14

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.

Mapa siedmiu 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

Cztery kolory wielokątów

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.

Whuber
źródło
Nie jestem pewien, czy to odpowiada na pytanie lub jakie jest pytanie, ale mimo wszystko jest to dobra odpowiedź.
nagytech
@Geoist Czy jest jakiś sposób, aby wyjaśnić ilustrację lub lepiej wyjaśnić problem?
whuber
6

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
ndimhypervol
źródło
2
Aby połączyć to z moim rozwiązaniem: twój countOverlapsodpowiada dwóm wierszom nbrhoods <- sapply(vertices, neighbors); degrees <- sapply(nbrhoods, length)w moim kodzie: degreesjest liczbą nakładania się. Oczywiście twój kod jest dłuższy, ponieważ odzwierciedla większość analiz GIS, które są brane za pewnik w moim rozwiązaniu: mianowicie, że najpierw identyfikujesz, które wielokąty zachodzą na siebie, i że na końcu używasz rozwiązania do generowania zestawów danych wieloboków. Dobrym pomysłem byłoby zamknięcie obliczeń teoretycznych na wykresach, więc jeśli kiedykolwiek znajdziesz lepszy algorytm kolorowania, łatwo byłoby go podłączyć.
whuber
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)
messamat
źródło
-3

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.

Alexandre Neto
źródło
2
ta metoda nie prowadzi do pożądanego produktu, który stanowi minimalną serię selekcji lub unikalnych klas funkcji oryginału, które się nie pokrywają. produkty zostaną wprowadzone do statystyki strefowej, dlatego kluczowe znaczenie ma zachowanie oryginalnej geometrii każdej cechy.
ndimhypervol
Masz rację, przepraszam. Nie zrozumiałem dobrze pytania. W takim przypadku, w zależności od wielkości rastra, normalnie przekształciłbym raster w tymczasową klasę elementów punktowych (każda komórka punkt) i wykonałbym przestrzenne połączenie między nim a warstwą wielokąta. Być może jest to bardzo uproszczone i nieprzyjazne pod względem wydajności podejście, ale działa, a nakładające się wielokąty nie będą stanowić problemu.
Alexandre Neto
Jeśli dobrze rozumiem, co rozumiesz przez to połączenie przestrzenne, twoje drugie rozwiązanie nadal nie zadziała, Alexandre, ponieważ istnieje relacja wiele do wielu między punktami a wielokątami. Niezależnie od tego, dla każdego dużego rastra takie podejście oparte na wektorze będzie wyjątkowo nieefektywne, a dla dużych rastrów niemożliwe.
whuber
@ whuber Masz rację, że jesteś bardzo powolnym procesem (Toke me około pół godziny z rastrem 4284 x 3009 i 2401 wielokątów, w dualcore 2.8Ghz, 3Gb RAM z perspektywą). Ale to działa, ponieważ już to przetestowałem. W sprzężeniu przestrzennym musisz użyć relacji jeden do jednego i agregować wartości rastrowe (jak średnia, Suma itp.). Rezultatem będzie warstwa wieloboku wektorowego podobna do oryginału, ale z nową kolumną z zagregowanymi wartościami rastrowymi, które przecinają każdy wielokąt. Nie jest to optymalne rozwiązanie, może to być przydatne dla osób o mniejszych umiejętnościach programistycznych (takich jak ja :-)).
Alexandre Neto,