¿Qué herramientas de generalización / suavizado de ráster están disponibles?

46

Tengo un DEM que me gustaría suavizar o generalizar para eliminar los extremos topográficos (cortar picos y valles de relleno). Idealmente, también me gustaría tener control sobre el radio o nivel de "borrosidad". Al final, necesitaré un conjunto de rásteres que van desde un poco borrosos a muy borrosos. (Teóricamente, el desenfoque sería una trama constante de la media aritmética de todos los valores).

¿Hay alguna herramienta o método que pueda usar (basado en Esri, GDAL, GRASS)? ¿Necesito hornear en casa mi propia rutina de desenfoque gaussiano ? ¿Podría usar un filtro de paso bajo (por ejemplo, el filtro de ArcGIS ) y, de ser así, ¿tendría que ejecutarlo varias veces para obtener el efecto de un radio grande?

Mike T
fuente
¿Qué pasa con solo exportar el ráster a un tamaño de celda más grande? ¿No resultaría esto también en un silenciamiento de extremos?
1
Sí, eso también reduciría los extremos (suponiendo que el remuestreo implícito implique algún tipo de promedio), pero es una forma terrible de suavizar un DEM: crearía una pequeña cantidad de bloques grandes. Por cierto, generalmente no es necesario exportar un ráster para hacer esto; La agregación y el remuestreo a un tamaño de celda diferente son operaciones básicas que generalmente se encuentran en el software basado en ráster.
whuber

Respuestas:

29

El desenfoque gaussiano es solo una media focal ponderada. Puede recrearlo con gran precisión con una secuencia de vecindad circular de corta distancia (sin ponderar) significa: esta es una aplicación del Teorema del límite central .

Tienes muchas opciones. El "filtro" es demasiado limitado, es solo para vecindarios de 3 x 3, así que no te molestes. La mejor opción para los DEM grandes es llevar el cálculo fuera de ArcGIS a un entorno que utilice transformaciones rápidas de Fourier: hacen los mismos cálculos focales pero (en comparación) lo hacen a una velocidad increíblemente rápida. (GRASS tiene un módulo FFT . Está destinado al procesamiento de imágenes, pero es posible que pueda ponerlo en servicio para su DEM si puede reescalarlo con una precisión razonable en el rango de 0..255). Salvo eso, al menos dos soluciones son vale la pena considerarlo:

  1. Cree un conjunto de pesos de vecindario para aproximar un desenfoque gaussiano para un vecindario considerable. Use pases sucesivos de este desenfoque para crear su secuencia de DEMs cada vez más suaves.

    (Los pesos se calculan como exp (-d ^ 2 / (2r)) donde d es la distancia (en celdas si lo desea) yr es el radio efectivo (también en celdas). Deben calcularse dentro de un círculo que se extiende al menos 3r . Después de hacerlo, divida cada peso por la suma de todos, de modo que al final sumen 1.)

  2. Alternativamente, olvide la ponderación; simplemente ejecute una media focal circular repetidamente. He hecho exactamente esto para estudiar cómo cambian las cuadrículas derivadas (como pendiente y aspecto) con la resolución de un DEM.

Ambos métodos funcionarán bien, y después de las primeras pasadas habrá poco para elegir entre los dos, pero hay rendimientos decrecientes: el radio efectivo de n medios focales sucesivos (todos con el mismo tamaño de vecindario) es solo (aproximadamente) el raíz cuadrada de n veces el radio de la media focal. Por lo tanto, para grandes cantidades de desenfoque, querrá comenzar de nuevo con un vecindario de gran radio. Si usa una media focal no ponderada, ejecute 5-6 pases sobre el DEM. Si usa pesos que son aproximadamente gaussianos, solo necesita una pasada: pero debe crear la matriz de peso.

Este enfoque tiene la media aritmética de la DEM como valor límite.

whuber
fuente
1
Si sus datos tienen picos, primero puede probar un filtro mediano ( en.wikipedia.org/wiki/Median_filter ) antes de aplicar un desenfoque más general como lo sugiere whuber.
MerseyViking
@ Mersey Esa es una excelente sugerencia. Nunca he visto un DEM con valores atípicos locales, pero tampoco he tenido que procesar un DEM sin procesar (como resultados LIDAR sin procesar). No puede hacer filtros de mediana con FFT, pero solo (generalmente) necesita un entorno de 3 x 3, por lo que es una operación rápida de todos modos.
Whuber
Gracias whuber. Debo admitir que solo he usado datos LiDAR preprocesados, pero hay algunos picos significativos en los datos SRTM que se beneficiarían de un filtro mediano. Sin embargo, tienden a tener 2 o 3 muestras de ancho, por lo que se necesitaría un filtro mediano más grande.
MerseyViking
@ Mersey Aún estás bien con un filtro mediano más grande de 5 x 5 o 7 x 7. Si estás contemplando (por ejemplo) un filtro de 101 x 101, ¡prepárate para esperar! También sugiere un punto importante que vale la pena elaborar: es una muy buena idea realizar un análisis exploratorio del DEM antes de hacer algo. Esto incluiría identificar picos (valores atípicos locales) y caracterizar sus tamaños y extensiones. ¡Desea asegurarse de que sean realmente artefactos (y no un fenómeno real) antes de eliminarlos con un filtro!
whuber
1
+1 para FFT en datos de elevación. De hecho, he hecho ese trabajo en el césped para datos NED de 32 bits para eliminar las rayas bidireccionales. Al final, esto también fue problemático porque reintrodujo el efecto de terrazas que plaga muchos otros DEM derivados del contorno.
Jay Guarneri
43

He estado explorando el enfoque signal.convolve de SciPy (basado en este libro de cocina ), y estoy teniendo un éxito realmente bueno con el siguiente fragmento:

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

Lo uso en otra función que lee / escribe geoatff float32 a través de GDAL (no es necesario volver a escalar a 0-255 bytes para el procesamiento de imágenes), y he estado usando tamaños de píxeles (por ejemplo, 2, 5, 20) y tiene salida realmente agradable (visualizada en ArcGIS con píxeles 1: 1 y rango mínimo / máximo constante):

Gaussian DTM

Nota: esta respuesta se actualizó para utilizar una función de procesamiento signal.fftconvolve basada en FFT mucho más rápida .

Mike T
fuente
1
+1 ¡Buena solución! No estoy seguro, pero es una buena apuesta que esa señal .volve utiliza FFT.
whuber
Estaba buscando un código borroso para una herramienta de costura automática que estoy escribiendo y me topé con esto. Buen trabajo @ MikeToews!
Ragi Yaser Burhum
@RagiYaserBurhum Me encantaría saber más sobre tu herramienta. MikeToews Gran respuesta y fragmento de código muy apreciado.
Jay Laura
@JayLaura Nada especial, solo escribo una herramienta para autostitch algunas imágenes que tomé con algunos amigos con un globo. Uso de las clases de Orfeo Toolbox orfeo-toolbox.org/SoftwareGuide/…
Ragi Yaser Burhum
2
@whuber al revisar esta rutina, no estaba usando FFT, pero ahora sí, y es mucho más rápido.
Mike T
4

Esto podría ser un comentario a la excelente respuesta de MikeT , si no fuera demasiado largo y complejo. He jugado mucho con él e hice un complemento QGIS llamado FFT Convolution Filters (todavía en etapa "experimental") basado en su función. Además de suavizar, el complemento también puede agudizar los bordes restando el ráster suavizado del original.

He mejorado un poco la función de Mike en el proceso:

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)

Las comprobaciones de validez son bastante evidentes, pero lo importante es lanzar y flotar. Antes de esto, la función hacía que las matrices de enteros fueran negras (solo ceros), debido a la división por la suma de los valores ( g / g.sum()).

Pavel V.
fuente
3

En QGIS, obtuve buenos resultados fácilmente al usar el filtrado de imágenes de Orfeo Toolbox . Es razonablemente rápido y el modo por lotes funciona bien. Se encuentran disponibles difusiones gaussianas, medias o anisotrópicas.

Tenga en cuenta que se Radiusrefiere al número de celdas, no a la distancia.

Aquí hay un ejemplo usando Smoothing (gaussian) :

  • Crudo:

    Sin filtro

  • Filtrado:

    filtrar

Tactopoda
fuente
1

Buena solución para el desenfoque gaussiano y la animación genial. Con respecto a la herramienta de filtro Esri mencionada anteriormente, esto es básicamente solo la herramienta Esri "Focal Statistics" codificada en un tamaño de 3x3. La herramienta de Estadísticas Focales le ofrece muchas más opciones sobre la forma de su filtro móvil, el tamaño y la estadística que desea ejecutar. http://desktop.arcgis.com/en/arcmap/latest/tools/spatial-analyst-toolbox/focal-statistics.htm

También puede hacer un filtro "irregular" donde pase su propio archivo de texto con pesos para usar en cada celda. El archivo de texto tiene tantas filas como desee en su área de filtro, con valores delimitados por espacios en blanco para las columnas. Supongo que siempre debe usar un número impar de filas y columnas, por lo que su celda objetivo está en el medio.

Creé una hoja de cálculo de Excel para jugar con diferentes pesos que acabo de copiar / pegar en este archivo. Debería lograr los mismos resultados que anteriormente si ajusta las fórmulas.

David A
fuente