46

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?

  • What about just exporting the raster to a larger cell size? Wouldn't this also result in a muting of extremes? – user14258 Jan 15 '13 at 16:18
  • 1
    Yes, that would also reduce extremes (assuming that the implicit resampling involves some form of averaging) but it's a terrible way to smooth a DEM: you would create a small number of large blocks. BTW, one usually does not need to export a raster to do this; aggregation as well as resampling to a different cellsize are basic operations usually found in raster-based software. – whuber Jan 15 '13 at 16:44
29

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.

  • 1
    If your data has spikes, you could try a median filter (en.wikipedia.org/wiki/Median_filter) first before applying a more general blur as suggested by whuber. – MerseyViking May 9 '11 at 8:44
  • @Mersey That's an excellent suggestion. I have never seen a DEM with local outliers, but then again I have never had to process a raw DEM (such as raw LIDAR results) either. You can't do median filters with FFT, but you only (usually) need a 3 x 3 neighborhood so it's a fast operation anyway. – whuber May 9 '11 at 14:49
  • Thanks whuber. I must admit I've only ever used pre-processed LiDAR data, but there are some significant spikes in SRTM data that would benefit from a median filter. They do tend to be 2 or 3 samples wide though, so a larger median filter would be needed. – MerseyViking May 9 '11 at 15:19
  • @Mersey You're still ok with a larger median filter of 5 x 5 or 7 x 7. If you're contemplating (say) a 101 x 101 filter, though, be prepared to wait! You also suggest an important point worth elaborating: it's a very good idea to perform an exploratory analysis of the DEM before doing anything. This would include identifying spikes (local outliers) and characterizing their sizes and extents. You want to be sure they're really artifacts (and not some real phenomenon) before you go about wiping them out with a filter! – whuber May 9 '11 at 17:12
  • 1
    +1 for FFT on elevation data. I've actually made that work in grass for 32bit NED data to remove bi-directional striping. In the end, this was also problematic because it re-introduced the terracing effect that plagues many other contour-derived DEMs. – Jay Guarneri May 1 '13 at 21:45
42

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.

  • 1
    +1 Nice solution! I don't know for sure, but it's a good bet that signal.convolve uses FFTs. – whuber Jun 14 '11 at 20:28
  • I was looking for some blurring code for an auto-stitching tool I am writing and stumbled upon this. Nice job @MikeToews! – Ragi Yaser Burhum Jan 15 '13 at 18:07
  • @RagiYaserBurhum Would love to hear more about your tool. MikeToews Great answer and much appreciated code snippet. – Jay Laura Jan 16 '13 at 2:55
  • @JayLaura Nothing special, just writing a tool to autostitch some images I took with some friends with a balloon. Using the Orfeo Toolbox classes orfeo-toolbox.org/SoftwareGuide/… – Ragi Yaser Burhum Jan 16 '13 at 5:37
  • 2
    @whuber upon revising this routine, it wasn't using FFT, but it is now, and is so much faster. – Mike T May 12 '14 at 1:28
4

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()).

3

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:

    No filter

  • Filtered:

    filter

1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.