Open a DEM, do a gaussian blur, then save it again Heavily based on: http://gis.stackexchange.com/questions/9431/what-raster-smoothing-generalization-tools-are-available

smoothDEM.py

#
# Open a DEM, do a gaussian blur, then save it again
# Heavily based on:
#   http://gis.stackexchange.com/questions/9431/what-raster-smoothing-generalization-tools-are-available
#
# Syntax for running snippet:
#   python smoothDem.py input.tif output.tif 5
#


from scipy import exp, mgrid, signal, row_stack, column_stack, tile, misc
from osgeo import gdal
import sys


SOURCE_FILE = sys.argv[1]
OUTPUT_FILE = sys.argv[2]
SIZE = int(sys.argv[3])


def gaussian_blur(in_array, size):

    # expand in_array to fit edge of kernel, by repeating edge values
    ny, nx = in_array.shape
    ary = row_stack((tile(in_array[0,:], size).reshape((size, nx)),
                     in_array,
                     tile(in_array[-1,:], size).reshape((size, nx))))
    arx = column_stack((tile(ary[:,0], size).reshape((size, ny + size*2)).T,
                        ary,
                        tile(ary[:,-1], size).reshape((size, ny + size*2)).T))
    
    # build kernel
    x, y = mgrid[-size:size+1, -size:size+1]
    g = exp(-(x**2/float(size) + y**2/float(size)))
    g = (g / g.sum()) #.astype(in_array.dtype)
    
    # do the gaussian blur over arx
    return signal.convolve(arx, g, mode='valid')


def createRasterArray():
    
    # Open source file
    source = gdal.Open(SOURCE_FILE)
    gt = source.GetGeoTransform()
    proj = source.GetProjection()
        
    # Get rasterarray
    band_array = source.GetRasterBand(1).ReadAsArray()
    blurredArray = gaussian_blur(band_array, SIZE)

    # Register driver
    driver = gdal.GetDriverByName('GTIFF')
    driver.Register()

    # Get source metadata
    cols = source.RasterXSize
    rows = source.RasterYSize
    bands = source.RasterCount
    band = source.GetRasterBand(1)
    datatype = band.DataType

    # Create output image
    output = driver.Create(OUTPUT_FILE, cols, rows, bands, datatype)
    output.SetGeoTransform(gt)
    output.SetProjection(proj)

    # Get band from newly created image
    outBand = output.GetRasterBand(1) 

    # Write to file
    outBand.WriteArray(blurredArray, 0, 0)
    print "Your file is saved!" 

    # Close all
    source = None # close raster
    output = None
    band = None
    outBand = None



createRasterArray()