[GIS] What raster smoothing/generalization tools are available

generalizationrastersmoothingsoftware-recommendations

I have a DEM that I would like to smooth or generalize to remove topographic extremes (chop off peaks and fill-in valleys). Ideally, I would also like to have control over the radius or level of "blurriness". In the end, I will need a set of rasters that range from slightly blurry to really blurry. (Theoretically, the blurriest would be a constant raster of the arithmetic mean of all values).

Are there any tools or methods I can use (based on Esri, GDAL, GRASS)?

Do I need to home bake my own Gaussian blur routine?

Could I use a low-pass filter (e.g. ArcGIS's filter), and if so, would I need to run it a bunch of times to get the effect of a large radius?

Best Answer

I've been exploring SciPy's signal.convolve approach (based on this cookbook), and am having some really nice success with the following snippet:

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

I use this in another function which reads/writes float32 GeoTIFFs via GDAL (no need to rescale to 0-255 byte for image processing), and I've been using attempting pixel sizes (e.g., 2, 5, 20) and it has really nice output (visualized in ArcGIS with 1:1 pixel and constant min/max range):

Gaussian DTM

Note: this answer was updated to use a much faster FFT-based signal.fftconvolve processing function.

Related Question