scipy tutorial question

[Xristos Xristoou]

hello

i want to follow this scipy tutorial http://www2.geog.ucl.ac.uk/~plewis/geogg122-2011-12/dem2.html to calculate slope/aspect and more.

first question,this is a good way to calculate slope/aspect with scipy ?
can i find more easy way for this calculate ?

second question and most important,how can i use mine DEM (raster image) inside
the tutorial code to calculate slope/aspect from mine DEM (raster image)?
and how to export that calculation to new tiff images ?
i thing so easy but i am very new.
i have and gdal package in my python

the tutorial code :


import numpy as np

def gaussianFilter(sizex,sizey=None,scale=0.333):
    '''
    Generate and return a 2D Gaussian function
    of dimensions (sizex,sizey)

    If sizey is not set, it defaults to sizex
    A scale can be defined to widen the function (default = 0.333)
    '''
    sizey = sizey or sizex
    x, y = np.mgrid[-sizex:sizex+1, -sizey:sizey+1]
    g = np.exp(-scale*(x**2/float(sizex)+y**2/float(sizey)))
    return g/g.sum()

def randomDEM(nx,sizex,max,min,sizey=None,ny=None):
    '''
    Generate a pseudo-DEM from random correlated noise.

    Parameters:
        nx  :       number of pixels of the output
    sizex   :       number of pixels of the smoothing filter
    max     :       max height in DEM
    min     :       min height in DEM

    Options:
    sizey   :       defaults to sizex
            ny      :       defaults to nx

    '''
    from scipy import signal
    ny = ny or nx
    sizey = sizey or sizex
    dem1 = np.random.rand(nx+2*sizex,ny+2*sizey)
    demSmooth = signal.convolve(dem1,gaussianFilter(sizex,sizey),mode='valid')
    demSmooth = (demSmooth - demSmooth.min())/(demSmooth.max() - demSmooth.min())
    demSmooth = demSmooth*(max-min)+min
    return demSmooth

def plot(array,file=None):
    '''
    Produce a plot of the 2D array

    Options:
    file    :       specify an output file
    '''
    import matplotlib.pylab as plt
    plt.imshow(array,interpolation='nearest')
    if file == None:
        plt.show()
    else:
        plt.savefile(file)
smoothDEM = randomDEM(300,20,100.,0.,sizey=5)
#plot(smoothDEM)
def padding(dem,size=1):
    '''
    Apply a border of size to a spatial dataset

    Return the padded data with the original centred in the array
    '''
    out = np.zeros([dem.shape[0]+2*size,dem.shape[1]+2*size])
    out[:,:] = np.max(dem)+1
    out[size:-size,size:-size] = dem
    return out

padDEM = padding(smoothDEM)
#plot(padDEM)
def localMin(dem):
    '''
    We wish to return the location of the minimum grid value
    in a neighbourhood.

    We assume the array is 2D and defined (y,x)

    We return wx,wy which are the cell displacements in x and y directions.

    '''
    wy = np.zeros_like(dem).astype(int)
    wx = np.zeros_like(dem).astype(int)
    winx = np.ones([3,3])
    for i in xrange(3):
        winx[:,i] = i - 1
    winy = winx.transpose()
    demp = padding(dem,size=1)
    for y in np.arange(1,demp.shape[0]-1):
        for x in np.arange(1,demp.shape[1]-1):
            win = demp[y-1:y+2,x-1:x+2]
            ww = np.where(win == np.min(win))
            whereX = winx[ww][0]
            whereY = winy[ww][0]
            wy[y-1,x-1] = whereY
            wx[y-1,x-1] = whereX
    return wx,wy

(wx,wy) = localMin(smoothDEM)
#plot(wx)


def grad2d(dem):
    '''
    Calculate the slope and gradient of a DEM
    '''
    from scipy import signal
    f0 = gaussianFilter(3)
    I = signal.convolve(dem,f0,mode='valid')
    f1 = np.array([[-1,0,1],[-2,0,2],[-1,0,1]])
    f2 = f1.transpose()
    g1 = signal.convolve(I,f1,mode='valid')
    g2 = signal.convolve(I,f2,mode='valid')
    slope = np.sqrt(g1**2 + g2**2)
    aspect = np.arctan2(g2,g1)
    return slope, aspect
slope,  aspect = grad2d(smoothDEM)
plot(slope)
plot(aspect)