What raster smoothing/generalization tools are available?

Question

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?

Answer 1

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

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

Answer 2

Gaussian blur is just a weighted focal mean. You can recreate it to high accuracy with a sequence of short-distance circular neighborhood (unweighted) means: this is an application of the Central Limit Theorem.

You have a lot of choices. "Filter" is too limited--it's only for 3 x 3 neighborhoods--so don't bother with it. The best option for large DEMs is to take the calculation outside of ArcGIS into an environment that uses Fast Fourier Transforms: they do the same focal calculations but (in comparison) they do it blazingly fast. (GRASS has an FFT module. It's intended for image processing but you might be able to press it into service for your DEM if you can rescale it with reasonable precision into the 0..255 range.) Barring that, two solutions at least are worth considering:

1. Create a set of neighborhood weights to approximate a Gaussian blur for a sizable neighborhood. Use successive passes of this blur to create your sequence of ever smoother DEMs.

   (The weights are computed as exp(-d^2/(2r)) where d is the distance (in cells if you like) and r is the effective radius (also in cells). They have to be computed within a circle extending out to at least 3r. After doing so, divide each weight by the sum of them all so at the end they sum to 1.)

2. Alternatively, forget the weighting; just run a circular focal mean repeatedly. I have done exactly this for studying how derived grids (like slope and aspect) change with the resolution of a DEM.

Both methods will work well, and after the first few passes there will be little to choose between the two, but there are diminishing returns: the effective radius of n successive focal means (all using the same neighborhood size) is only (approximately) the square root of n times the radius of the focal mean. Thus, for huge amounts of blurring, you will want to begin over again with a large-radius neighborhood. If you use an unweighted focal mean, run 5-6 passes over the DEM. If you use weights that are approximately Gaussian, you need only one pass: but you have to create the weight matrix.

This approach indeed has the arithmetic mean of the DEM as a limiting value.

Answer 3

This could be a comment to MikeT's excellent answer, if it wasn't too long and too complex. I've played with it a lot and made a QGIS plugin named FFT Convolution Filters (in "experimental" stage yet) based on his function. Besides smoothing, the plugin can also sharpen edges by subtracting the smoothed raster from the original one.

I've upgraded Mike's function a little in the process:

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)

The validity checks are quite self-evident, but what's important is casting to float and back. Before this, the function made integer arrays black (zeros only), because of the dividing by the sum of the values (g / g.sum()).

Answer 4

In QGIS, I got good results easily by using Orfeo Toolbox Image filtering. It's reasonable fast and batch mode works fine. Gaussian, mean, or anisotropic diffusions are available.

Note that Radius refers to number of cells, not distance.

Here is an example using Smoothing(gaussian):

• Raw:

  

• Filtered:

  

Answer 5

Nice solution for the Gaussian blur and cool animation. Regarding the Esri Filter tool mentioned above, that is basically just the Esri "Focal Statistics" tool hard-coded to a 3x3 size. The Focal Statistics tool gives you a lot more options on the shape of your moving filter, the size, and the statistic you want to run. http://desktop.arcgis.com/en/arcmap/latest/tools/spatial-analyst-toolbox/focal-statistics.htm

You can also make an "irregular" filter where you pass in your own text file with weights to use for each cell. The text file has as many rows as you want in your filter area, with whitespace-delimited values for the columns. I guess you should always use an odd numbers of rows and columns, so your target cell is in the middle.

I created an excel spreadsheet to play with different weights that I just copy/paste into this file. It should achieve the same results as above if you adjust the formulas.