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)
--
https://mail.python.org/mailman/listinfo/python-list
