Contando valores de pixels consecutivos para um conjunto de rasters usando o ArcGIS Spatial Analyst?

23

Estou usando o ArcGIS 10 com o Spatial Analyst.

Eu tenho um conjunto de rasters (8 no total) que apenas todos contêm 1 ou 0 para cada célula. Cada varredura representa um ano diferente de dados. Para argumentos, ano 1 a 8.

Eu posso adicionar todos os rasters juntos, o que me dará uma grade final com valores variando de 0 a 8. Um 8 indicando que a célula era constantemente 1 para o conjunto de rasters (todos os anos).

Gostaria de descobrir para cada célula o maior número consecutivo de 1's.

Assim, por exemplo, a grade total pode registrar para uma única célula um valor de digamos 5, mas nas 8 grades essa célula possui o maior número consecutivo de 1's igual a 3. Ou outra maneira de expressar isso é por 3 anos que a célula era 1 então começou a oscilar entre zeros e uns.

Minhas habilidades de processamento de varredura não são tão boas quanto minha habilidade de processamento de vetores e eu dei uma boa olhada no arquivo de ajuda da ESRI, mas não consigo descobrir como alguém conseguiria isso usando as ferramentas de processamento geográfico prontas para uso?

Alguma ideia?

Hornbydd
fonte
1
Essa é realmente uma análise bem legal. Como sempre, há mais de uma maneira de fazer o que você está tentando fazer. Eu acho que você precisará fazer alguma programação para percorrer todas as combinações.
MLowry
1
Um comentário geral (inspirado por esta observação do @MLowry): faça um voto positivo sempre que parecer interessante ou claramente articulado. Boas perguntas direcionam tudo em nosso site; por favor, façamos o possível para recompensar quem pedir!
whuber

Respostas:

14

Como esta é uma operação local, vamos descobrir como fazer isso para uma única célula: o Map Algebra cuidará do resto.

Primeira nota que a ordem dos rasters obviamente importa. Portanto, as estatísticas de célula de disparo único, como uma soma de células, não o farão.

Se encontrarmos uma sequência como 01110101 em uma determinada célula, processaremos isso do início ao fim e

  1. Comece com uma contagem em zero.

  2. Aumente a contagem sempre que encontrarmos 1.

  3. Redefina a contagem sempre que encontrarmos um 0, depois de salvar a última contagem .

  4. No final, faça a contagem máxima salva (incluindo a contagem final).

A etapa 1 é implementada com uma grade zero constante. Os passos 2 e 3 dependem do que encontramos: é, portanto, uma operação condicional . O passo 4 é claramente um máximo local. Codificaríamos isso, então, um pouco mais formalmente como:

count = 0
result = 0
For each value:
    If (value==1):
        count=count+1
    else
        result = max(result, count)
        count=0
result = max(result, count)

É melhor fazer isso com um script Python quando você tem muitas grades, mas com oito não é oneroso desenrolar o loop e escrever as etapas manualmente. Isso revela um pequeno problema: result=max(longest,count)é um tipo de "efeito colateral" difícil de codificar com operações de varredura. (Mas isso pode ser feito, conforme mostrado na segunda solução abaixo.) Também é ineficiente, porque adiciona um cálculo extra a cada etapa. Portanto, modificamos um pouco a abordagem, com o objetivo de adiar a maxoperação até o final. Isso exigirá salvar uma contagem separada em cada estágio.

Ao passar por esse processo, também encontrei um atalho para o primeiro passo. Isso leva à seguinte solução, que embora um pouco longa e com muita memória RAM, é simples e envolve etapas executadas rapidamente:

result1 = "grid1"
result2 = con("grid2"==1, "result1"+1, 0)
result3 = con("grid3"==1, "result2"+1, 0)
result4 = con("grid4"==1, "result3"+1, 0)
result5 = con("grid5"==1, "result4"+1, 0)
result6 = con("grid6"==1, "result5"+1, 0)
result7 = con("grid7"==1, "result6"+1, 0)
result8 = con("grid8"==1, "result7"+1, 0)
CellStatistics(["result1", "result2", "result3", "result4", "result5", "result6", "result7" "result8"], "max")

A sintaxe real varia de acordo com a sua versão do ArcMap. (Por exemplo, CellStatisticsé novo na versão 10, acredito, mas uma operação máxima local sempre esteve disponível.)

No exemplo com a entrada 01110101, a sequência de grades "result *" conteria os valores 0, 1, 2, 3, 0, 1, 0, 1, portanto, no final CellStatistics, retornaria 3, o comprimento da cadeia mais longa de 1's.

Se a RAM for escassa, a solução poderá ser modificada para reutilizar os resultados intermediários, a um custo de aproximadamente o dobro do tempo de execução:

result = "grid1"
temp = con("grid2"==1, "result"+1, 0)
result = CellStatistics[["temp", "result"], "max"]
temp = con("grid3"==1, "temp"+1, 0)
result = CellStatistics[["temp", "result"], "max"]
...
temp = con("grid8"==1, "temp"+1, 0)
CellStatistics[["temp", "result"], "max"]

No exemplo com a entrada 01110101, os valores ("temp", "result") seriam (NoData, 0) após a primeira linha e após cada par de operações ("Con", "CellStatistics") os valores seriam (1 , 1), (2, 2), (3, 3), (0, 3), (1, 3), (0, 3), (1, 3). Mais uma vez, o valor final é 3.

O padrão regular de expressões de Álgebra de Mapa em qualquer solução indica como codificar o algoritmo como um loop em um script, alterando os índices conforme apropriado a cada iteração.

whuber
fonte
Pode haver um erro de digitação no bloco de código do tipo: count = count = 1 provavelmente deve ser count = count + 1
MLowry
1
@ ML Obrigado (bons olhos!): Está corrigido agora. É difícil tornar o pseudocódigo absolutamente correto; A revisão humana é um trunfo real na busca de erros. Além disso, embora eu não tenha testado as soluções no ArcGIS, implementei a primeira solução em R, por isso tenho certeza de que a abordagem está correta.
whuber
1
"whuber" mais uma vez você é o homem que sabe! Deus ajude o resto de nós se você for atropelado por um ônibus! Sua abordagem inicial do Python foi a direção que eu estava pensando, mas eu sei que com os rasters, muitas vezes é possível fazer tudo com muito mais clareza que você provou. Se você se encontrar na Grã-Bretanha, seria uma honra comprar uma caneca da nossa melhor cerveja plana em temperatura ambiente! :)
Hornbydd
Obrigado, Duncan: mas confira a excelente solução de Andy Harfoot!
whuber
14

Apenas conversando sobre isso e imaginando se você poderia abordar o problema tratando as grades de entrada como um fluxo binário. Isso permitiria combiná-los para fornecer um número inteiro de resumo exclusivo para a sequência - ou seja, 01110101 = 117. Esse valor pode ser reclassificado para fornecer o número máximo de 1s consecutivos.

Aqui está um exemplo mostrando uma maneira de combinar oito grades:

2*(2*(2*(2*(2*(2*(2*"g8" + "g7") + "g6") + "g5") + "g4") + "g3") + "g2") + "g1"

As operações bit a bit também podem ser pressionadas em serviço para esta etapa. Como alternativa, você pode usar a combinação seguida por um cálculo de campo. (O cálculo do campo terá uma expressão semelhante à anterior.)

A tabela de reclassificação deve fornecer comprimentos máximos de execução para todos os valores entre 00000000B = 0 e 11111111B = 255. Em ordem, aqui estão eles:

0, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 1, 1, 1, 2, 1, 1, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 5, 6, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 1, 1, 1, 2, 1, 1, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 3, 4, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 6, 7, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 1, 1, 1, 2, 1, 1, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 5, 6, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 3, 4, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 7, 8

Essa abordagem é limitada a cerca de 20 grades no ArcGIS: usar mais do que isso pode criar uma tabela de atributos complicada. ( Combineestá especificamente limitado a 20 grades).

Andy Harfoot
fonte
Um +1 enorme: essa é uma ideia muito boa . (A única limitação é que, quando mais de 31 grades estiverem envolvidas, você ficará sem bits para usá-las.) Tomei a liberdade de aprofundar um pouco sua ideia para que outros possam ver como é fácil implementar.
whuber
3

Você já pensou em alterar os valores de 0 e 1 para valores com a potência de 2 (1,2,4,8,16,32). Ao combinar as 8 grades, você obterá valores exclusivos para cada célula, que fornecerão informações consecutivas (ou seja: um valor de 3 significa os anos 1 e 2, onde um valor de 54 seria os anos de 6 a 8).

Apenas um pensamento

Ryan Garnett
fonte
Foi exatamente isso que @Andy Harfoot sugeriu algumas horas antes, Ryan. :-)
whuber
Obrigado e desculpe. Eu li isso no meu telefone de férias.
Ryan Garnett