Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question

4 Answers 4

up vote 17 down vote accepted

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.

share|improve this answer
    
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

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.

share|improve this answer
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
1  
@whuber upon revising this routine, it wasn't using FFT, but it is now, and is so much faster. –  Mike T May 12 at 1:28

I had success in smoothing a hillshade (reducing the pixeled effect) for map production using GDAL Warp method mentioned in this post

Ando

share|improve this answer
    
Normally smoothing is an unwanted side-effect of warping and warping algorithms are designed to minimize it. I suspect your hillshade appeared "smoother" because the warping may have resampled the grid to a smaller (relative) cellsize. In this particular application (of making a hillshade look nicer) that's a clever idea, but I don't think it generalizes to smoothing DEMs or grids in general. –  whuber May 9 '11 at 15:05
    
@whuber Yes the hillshade appeared smoother because I altered the cell size to a smaller value. Sorry I didn't read the question clearly enough to see that Mike wants to reduce/smooth the extreme values in the dataset!!! Whooops!!! –  Ando May 10 '11 at 22:19

What about just exporting the raster to a larger cell size? Wouldn't this also result in a muting of extremes?

share|improve this answer
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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