Quais ferramentas de suavização / generalização de varredura estão disponíveis?

46

Eu tenho um DEM que gostaria de suavizar ou generalizar para remover extremos topográficos (cortar picos e vales de preenchimento). Idealmente, eu também gostaria de ter controle sobre o raio ou o nível de "embaçamento". No final, precisarei de um conjunto de rasters que variam de levemente desfocado a realmente desfocado. (Teoricamente, o desfoque seria uma varredura constante da média aritmética de todos os valores).

Existem ferramentas ou métodos que posso usar (com base em Esri, GDAL, GRASS)? Preciso fazer em casa minha própria rotina de desfoque gaussiano ? Eu poderia usar um filtro passa-baixo (por exemplo, filtro do ArcGIS ) e, se sim, precisaria executá-lo várias vezes para obter o efeito de um raio grande?

Mike T
fonte
Que tal exportar a varredura para um tamanho de célula maior? Isso também não resultaria em silenciamento de extremos?
1
Sim, isso também reduziria extremos (assumindo que a reamostragem implícita envolva alguma forma de média), mas é uma maneira terrível de suavizar um DEM: você criaria um pequeno número de blocos grandes. BTW, normalmente não é necessário exportar uma varredura para fazer isso; agregação e reamostragem para um tamanho de célula diferente são operações básicas geralmente encontradas em software baseado em varredura.
whuber

Respostas:

29

O desfoque gaussiano é apenas uma média focal ponderada. Você pode recriá-lo com alta precisão com uma sequência de meios circulares de vizinhança de curta distância (sem peso): esta é uma aplicação do Teorema do Limite Central .

Você tem muitas opções. O "filtro" é muito limitado - é apenas para bairros 3 x 3 - por isso não se preocupe. A melhor opção para os grandes DEMs é levar o cálculo para fora do ArcGIS para um ambiente que use as transformações rápidas de Fourier: eles fazem os mesmos cálculos focais, mas (em comparação) eles fazem isso de uma maneira incrivelmente rápida. (O GRASS possui um módulo FFT . Destina-se ao processamento de imagens, mas você pode colocá-lo em serviço no seu DEM se for possível redimensioná-lo com razoável precisão na faixa de 0..255.) Exceto isso, pelo menos duas soluções Vale a pena considerar:

  1. Crie um conjunto de pesos de vizinhança para aproximar um desfoque gaussiano de uma vizinhança considerável. Use passes sucessivos desse desfoque para criar sua sequência de DEMs cada vez mais suaves.

    (Os pesos são calculados como exp (-d ^ 2 / (2r)) onde d é a distância (nas células, se desejar) e r é o raio efetivo (também nas células). Eles devem ser calculados dentro de um círculo que se estende até pelo menos 3. Depois de fazer isso, divida cada peso pela soma de todos, de modo que no final eles somam 1.)

  2. Como alternativa, esqueça a ponderação; basta executar uma média focal circular repetidamente. Fiz exatamente isso estudando como as grades derivadas (como inclinação e aspecto) mudam com a resolução de um DEM.

Ambos os métodos funcionarão bem e, após os primeiros passes, haverá pouco a escolher entre os dois, mas há retornos decrescentes: o raio efetivo de n meios focais sucessivos (todos usando o mesmo tamanho de vizinhança) é apenas (aproximadamente) o raiz quadrada de n vezes o raio da média focal. Assim, para grandes quantidades de desfoque, você desejará começar de novo com uma vizinhança de grande raio. Se você usar uma média focal não ponderada, execute 5-6 passagens sobre o DEM. Se você usar pesos aproximadamente gaussianos, precisará de apenas uma passagem: mas precisará criar a matriz de pesos.

Essa abordagem realmente tem a média aritmética do DEM como um valor limitador.

whuber
fonte
1
Se seus dados apresentarem picos, tente um filtro mediano ( pt.wikipedia.org/wiki/Median_filter ) antes de aplicar um desfoque mais geral, conforme sugerido pelo whuber.
MerseyViking
@ Jersey Essa é uma excelente sugestão. Eu nunca vi um DEM com discrepantes locais, mas também nunca tive que processar um DEM bruto (como resultados brutos do LIDAR) também. Você não pode fazer filtros medianos com a FFT, mas apenas (geralmente) precisa de uma vizinhança 3 x 3, portanto, é uma operação rápida de qualquer maneira.
whuber
Obrigado whuber. Devo admitir que só usei dados LiDAR pré-processados, mas há alguns picos significativos nos dados SRTM que se beneficiariam de um filtro mediano. Eles costumam ter 2 ou 3 amostras de largura, portanto, seria necessário um filtro mediano maior.
MerseyViking
@ Jersey Você ainda está bem com um filtro mediano maior de 5 x 5 ou 7 x 7. Se você estiver contemplando (digamos) um filtro 101 x 101, esteja preparado para esperar! Você também sugere um ponto importante que vale a pena elaborar: é uma boa idéia realizar uma análise exploratória do DEM antes de fazer qualquer coisa. Isso incluiria a identificação de picos (outliers locais) e a caracterização de seus tamanhos e extensões. Você quer ter certeza de que eles são realmente artefatos (e não um fenômeno real) antes de acabar com eles com um filtro!
whuber
1
+1 para FFT em dados de elevação. Na verdade, eu fiz esse trabalho na grama para dados NED de 32 bits para remover as faixas bidirecionais. No final, isso também foi problemático porque reintroduziu o efeito de aterramento que afeta muitos outros DEMs derivados de contorno.
Jay Guarneri
43

Eu estive explorando a abordagem signal.convolve do SciPy (com base neste livro de receitas ) e estou obtendo um sucesso muito bom com o seguinte trecho:

import numpy as np
from scipy.signal import fftconvolve

def gaussian_blur(in_array, size):
    # expand in_array to fit edge of kernel
    padded_array = np.pad(in_array, size, 'symmetric')
    # build kernel
    x, y = np.mgrid[-size:size + 1, -size:size + 1]
    g = np.exp(-(x**2 / float(size) + y**2 / float(size)))
    g = (g / g.sum()).astype(in_array.dtype)
    # do the Gaussian blur
    return fftconvolve(padded_array, g, mode='valid')

Eu uso isso em outra função que lê / grava GeoTIFFs float32 via GDAL (não é necessário redimensionar para 0-255 bytes para processamento de imagem), e eu tenho tentado tamanhos de pixel (por exemplo, 2, 5, 20) e ele tem saída realmente agradável (visualizada no ArcGIS com 1: 1 pixel e alcance mínimo / máximo constante):

Gaussian DTM

Nota: esta resposta foi atualizada para usar uma função de processamento signal.fftconvolve baseada em FFT muito mais rápida .

Mike T
fonte
1
+1 boa solução! Não sei ao certo, mas é uma boa aposta que o signal.convolve use FFTs.
whuber
Eu estava procurando algum código de desfoque para uma ferramenta de costura automática que estou escrevendo e me deparei com isso. Bom trabalho @MikeToews!
Ragi Yaser Burhum
@RagiYaserBurhum Gostaria de saber mais sobre sua ferramenta. MikeToews Ótima resposta e snippet de código muito apreciado.
Jay Laura
@ JayLaura Nada de especial, basta escrever uma ferramenta para fixar automaticamente algumas imagens que tirei com alguns amigos com um balão. Usando as classes Orfeo Toolbox orfeo-toolbox.org/SoftwareGuide/...
Ragi Yaser Burhum
2
@whuber ao revisar essa rotina, ela não estava usando o FFT, mas agora é e é muito mais rápido.
Mike T
4

Isso poderia ser um comentário à excelente resposta de MikeT , se não fosse muito longa e muito complexa. Eu brinquei muito com ele e criei um plugin QGIS chamado FFT Convolution Filters (no estágio "experimental" ainda) com base em sua função. Além de suavizar, o plugin também pode afiar as bordas subtraindo a varredura suavizada da original.

Atualizei a função de Mike um pouco no processo:

def __gaussian_blur1d(self, in_array, size):
        #check validity
        try:
            if 0 in in_array.shape:
                raise Exception("Null array can't be processed!")
        except TypeError:
            raise Exception("Null array can't be processed!")
        # expand in_array to fit edge of kernel
        padded_array = np.pad(in_array, size, 'symmetric').astype(float)
        # build kernel
        x, y = np.mgrid[-size:size + 1, -size:size + 1]
        g = np.exp(-(x**2 / float(size) + y**2 / float(size)))
        g = (g / g.sum()).astype(float)
        # do the Gaussian blur
        out_array = fftconvolve(padded_array, g, mode='valid')
        return out_array.astype(in_array.dtype)

As verificações de validade são bastante evidentes, mas o importante é lançar para flutuar e voltar. Antes disso, a função tornava as matrizes inteiras pretas (somente zeros), devido à divisão pela soma dos valores ( g / g.sum()).

Pavel V.
fonte
3

No QGIS, obtive bons resultados com facilidade usando o filtro de imagem do Orfeo Toolbox . É rápido e o modo em lote é razoável. Difusões gaussianas, médias ou anisotrópicas estão disponíveis.

Observe que Radiusse refere ao número de células, não à distância.

Aqui está um exemplo usando Suavização (gaussiana) :

  • Cru:

    Sem filtro

  • Filtrado:

    filtro

Tactopoda
fonte
1

Ótima solução para o borrão gaussiano e animação interessante. Com relação à ferramenta Esri Filter mencionada acima, essa é basicamente a ferramenta Esri "Estatísticas Focais" codificada para um tamanho 3x3. A ferramenta Estatísticas focais oferece muito mais opções sobre a forma do filtro em movimento, o tamanho e a estatística que você deseja executar. http://desktop.arcgis.com/en/arcmap/latest/tools/spatial-analyst-toolbox/focal-statistics.htm

Você também pode criar um filtro "irregular" onde passa seu próprio arquivo de texto com pesos para usar em cada célula. O arquivo de texto possui quantas linhas você deseja na sua área de filtro, com valores delimitados por espaço em branco para as colunas. Eu acho que você deve sempre usar um número ímpar de linhas e colunas, para que sua célula de destino esteja no meio.

Criei uma planilha do Excel para jogar com pesos diferentes que acabei de copiar / colar neste arquivo. Deverá alcançar os mesmos resultados acima, se você ajustar as fórmulas.

David A
fonte