Explodindo sobreposição a novos polígonos não sobrepostos?

10

Dado vários polígonos que se sobrepõem de várias maneiras, eu gostaria de exportar desses recursos todos os polígonos que não se sobrepõem aos outros, iterativamente.

O produto seria um número de recursos sem sobreposição que, quando somados, compõem o original.

Os produtos poderiam então ser usados ​​como entrada para as Estatísticas Zonais, e isso seria muito mais rápido do que repetir as Estatísticas Zonais sobre cada polígono.

Eu tenho tentado codificar isso no ArcPy sem sucesso.

O código para fazer isso já existe?

ndimhypervol
fonte
Você quer dizer que deseja 'achatar' os dados em um conjunto topologicamente correto?
Nagytech 26/08/12
O @Geoist ZonalStats requer polígonos que não se sobrepõem. Quando você tem uma coleção sobreposta, a solução óbvia, mas ineficiente, é fazer um loop sobre as políticas e calcular as estatísticas das zonas uma a uma. Seria mais eficiente selecionar um subconjunto de políticas não sobrepostas, aplicar estatísticas zonais a elas e iterar. A pergunta pergunta como fazer essas seleções com eficiência.
whuber
whuber - Eu acho que o @Geoist está sugerindo a criação de um conjunto de polígonos não sobrepostos a partir das interseções dos polígonos de entrada. Veja esta imagem - (não é possível postar imagens nos comentários?). A entrada está à esquerda. Toda a região é coberta por três polígonos, cada um dos quais cruza os outros. Os únicos subconjuntos que não se sobrepõem são os singletons e estes não atendem aos requisitos de Gotanuki de que a união preencha o espaço. Acho Geoist está sugerindo a criação do conjunto de regiões não-interseção à direita, que é válido para zonalstats
Llaves
Eu acho que há alguma confusão sobre qual deve ser o produto final. Você poderia dar um exemplo? Minha interpretação é que você gostaria que a saída fosse uma seleção de polígonos que não se sobrepõem - enquanto descarta ou dissolve os polígonos restantes. Você está trabalhando com uma ou várias classes de recursos?
Aaron
11
Parece-me que @gotanuki está querendo criar o número mínimo de classes de recurso que contêm apenas não sobreposição polígonos de uma classe de recurso polígono com polígonos sobrepostos
PolyGeo

Respostas:

14

Este é um problema de coloração gráfica .

Lembre-se de que a coloração de um gráfico é uma atribuição de uma cor aos vértices de um gráfico, de modo que nenhum dois vértices que compartilham uma aresta também tenham a mesma cor. Especificamente, os vértices (abstratos) do gráfico são os polígonos. Dois vértices são conectados a uma aresta (não direcionada) sempre que se cruzam (como polígonos). Se tomarmos alguma solução para o problema - que é uma sequência de (digamos, k ) coleções separadas dos polígonos - e atribuirmos uma cor única a cada coleção da sequência, obteremos uma coloração k do gráfico . É desejável encontrar um pequeno k .

Esse problema é bastante difícil e permanece sem solução para gráficos arbitrários. Considere uma solução aproximada que seja simples de codificar. Um algoritmo seqüencial deve fazer. O algoritmo de Galês-Powell é uma solução gananciosa baseada em uma ordem decrescente dos vértices por grau. Traduzido para o idioma dos polígonos originais, primeiro classifique os polígonos em ordem decrescente do número de outros polígonos sobrepostos. Trabalhando em ordem, dê ao primeiro polígono uma cor inicial. Em cada etapa sucessiva, tente colorir o próximo polígono com uma cor existente: ou seja, escolha uma cor que não sejajá usado por qualquer um dos vizinhos desse polígono. (Existem várias maneiras de escolher entre as cores disponíveis; tente a que foi menos usada ou escolha uma aleatoriamente.) Se o próximo polígono não puder ser colorido com uma cor existente, crie uma nova cor e pinte-a com ela.

Depois de obter uma coloração com um pequeno número de cores, execute as estatísticas zonais cor por cor: por construção, você garante que nenhum polígono de uma determinada cor se sobreponha.


Aqui está um exemplo de código R. (O código Python não seria muito diferente.) Primeiro, descrevemos sobreposições entre os sete polígonos mostrados.

Mapa de sete polígonos

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)

Ou seja, os polígonos 1 e 2 se sobrepõem, e os polígonos 2 e 3, 3 e 4, ..., 1 e 7.

Classifique os vértices por grau 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))]

Um algoritmo de coloração sequencial (bruto) usa a cor disponível mais antiga ainda não usada por qualquer polígono sobreposto:

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
}

Inicialize as estruturas de dados ( colorse color.next) e aplique o algoritmo:

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

Divida os polígonos em grupos de acordo com a cor:

split(vertices, colors)

A saída neste exemplo usa quatro cores:

$`1`
[1] 2 4

$`2`
[1] 3 6 7

$`3`
[1] 5

$`4`
[1] 1

Quatro cores dos polígonos

Ele particionou os polígonos em quatro grupos não sobrepostos. Nesse caso, a solução não é ótima ({{3,6,5}, {2,4}, {1,7}} são três cores para este gráfico). Em geral, a solução que obtém não deve ser tão ruim assim.

whuber
fonte
Não tenho certeza se isso responde à pergunta, ou qual é a pergunta, mas é uma boa resposta, no entanto.
Nagytech 28/08/12
@ Geoist Existe alguma maneira de tornar a ilustração mais clara ou explicar melhor o problema?
whuber
6

A metodologia recomendada pelo #whuber me inspirou a tomar uma nova direção, e aqui está minha solução arqueada, em duas funções. O primeiro, chamado countOverlaps, cria dois campos, "overlaps" e "ovlpCount", para gravar para cada poli que polys se sobrepõem a ele e quantas sobreposições ocorreram. A segunda função, explodeOverlaps, cria um terceiro campo, "expl", que fornece um número inteiro exclusivo para cada grupo de polys não sobrepostas. O usuário pode exportar novos fc com base nesse campo. O processo está dividido em duas funções, porque acho que a ferramenta countOverlaps pode ser útil por si só. Por favor, desculpe a negligência do código (e a convenção de nomenclatura descuidada), pois é bastante preliminar, mas funciona. Verifique também se o "idName" campo representa um campo com IDs exclusivos (testados apenas com IDs inteiros). Obrigado whuber por me fornecer a estrutura necessária 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
ndimhypervol
fonte
2
Para conectar isso à minha solução: o seu countOverlapscorresponde às duas linhas nbrhoods <- sapply(vertices, neighbors); degrees <- sapply(nbrhoods, length)do meu código: degreesé a contagem de sobreposições. É claro que seu código é mais longo porque reflete a maior parte da análise GIS que é dada como certa em minha solução: a saber, que você primeiro identifica quais polígonos se sobrepõem e que, no final, usa a solução para gerar conjuntos de dados de polígonos. Seria uma boa idéia encapsular os cálculos teóricos dos grafos; portanto, se você encontrar um algoritmo de coloração melhor, será fácil
inseri
1

Já faz um tempo, mas eu usei esse código para meu próprio aplicativo e ele está funcionando muito bem - obrigado. Reescrevi parte dele para atualizá-lo, aplicá-lo às linhas (com tolerância) e acelerar significativamente (abaixo - estou executando-o em 50 milhões de buffers que se cruzam e leva apenas algumas 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)
messamat
fonte
-3

Nesse caso, geralmente uso o seguinte método:

  • Passe a classe Feature através de uma UNION; (Ele quebra os polígonos em todas as suas interseções)
  • Adicione os campos X, Y e Área e calcule-os;
  • Dissolva o resultado pelos campos X, Y, Área.

Acredito que o resultado será o desejado e você pode até contar o número de sobreposições. Não sei se, em termos de desempenho, será melhor para você ou não.

Alexandre Neto
fonte
2
esse método não leva você ao produto desejado, que é uma série mínima de seleções ou classes de recurso exclusivas do original que não se sobrepõem. os produtos serão alimentados em estatísticas zonais e, portanto, é vital manter a geometria original de cada recurso.
Ndimhypervol
Você está certo, desculpe. Não entendi bem a pergunta. Nesse caso, e dependendo do tamanho da varredura, eu normalmente converteria a varredura em uma classe de recurso de ponto temporário (cada célula por ponto) e executaria uma junção espacial entre ela e a camada de polígono. Talvez seja uma abordagem muito simplista e hostil ao desempenho, mas funcione e os polígonos sobrepostos não lhe causem problemas.
Alexandre Neto
Se eu entendi direito o que você quer dizer com essa junção espacial, sua segunda solução ainda não funcionará, Alexandre, porque há uma relação de muitos para muitos entre os pontos e os polígonos. Independentemente disso, para qualquer varredura considerável, essa abordagem baseada em vetor será extremamente ineficiente e, para rasters grandes, será impossível executar.
whuber
@whuber Você está certo quanto a ser um processo muito lento (fale comigo por meia hora com uma varredura de 4284 x 3009 e 2401 polígonos, em uma RAM dualcore de 2,8 Ghz e 3Gb com vista). Mas funciona, como eu já testei. Na Junção Espacial, você deve usar uma relação de um para um e agregar os valores de varredura (como média, Soma, etc ...). O resultado será uma camada de polígono vetorial semelhante ao original, mas com uma nova coluna com os valores de varredura agregados que cruzam cada polígono. Não sendo uma solução ideal, isso pode ser útil para alguém com menos habilidades de programação (como eu :-)).
Alexandre Neto