Gully Erosion Model
From GeoMod
Gully erosion code from Ekta Amar's Class Project that has been modified to allow user interaction.
- User must press button to get to next timestep.
- User can set initial slope.
To run you must have the following 3 files in the same directory (with the names given here). You run the gully.py file; the other two files are imported by this program.
Contents |
[edit]
Project documentation
For details regarding the software and the equations employed, see Ekta's Class Project.
- The pictorial description of project outline provides a good overview of what the model is trying to accomplish.
- The literature review.
- Description of Features of the Model, including the model parameters.
- Documentation of the functions in the program.
- Images showing simulation results including the drainage network and the spatial distribution of erosion.
[edit]
Main program (gully.py)
from visual import *
from random import *
from SAREA_6_tmp import *
from uControls import *
from visual.graph import *
import math
#WINDOW
##origx = 0
##origy = 0
w = 704 #+4+4
h = 576 #+24+4
scene.width=w
scene.height=h
scene.x = 0
scene.y = 0
# Define parameters for creating the hypothetical topography
nx = 25
lx = 25.0
ny = 25
ly = 25.0
dx = lx/nx # Width of each cell in x direction
dy = ly/ny # Width of each cell in y direction
# creating land area as per the parameters defined
sourceArea = area (nx, ny, dx, dy)
# Defining general slope to the topography
initial_slope = 0.5
sourceArea.slope(dydz=initial_slope)
#controls
cwin = display(title="Controls", width=300, height=300,
x=0, y=h+24)
tstep_but = uSwitch(title="Timestep")
raise_slope_but = uSwitch(title="Increase slope", pos=(4,2.0))
lower_slope_but = uSwitch(title="Lower slope", pos=(4,-2.0))
pick = None
tt = 0 # timestep
while 1:
if cwin.mouse.events:
#print "event"
m2 = cwin.mouse.getevent() # obtain drag or drop event
#DO TIME STEP
if m2.pick == tstep_but.button and m2.release:
tstep_but.switch()
tt = tt+1
print tt
tstep_but.val_label.text = str(tt)
#rate(1) # Specifying the slow down rate
sourceArea.update_topo() # To update the elevation of the topography
sourceArea.catchment(tt) # To compute the lowest elevation cell in a array of 3x3
sourceArea.watershed(tt) # To compute the watershed basin for all cells
sourceArea.localSlope() # To compute the local slope for all cells
sourceArea.sediTrans(tt) # To compute the sediment transport capacity of all cells
sourceArea.eroht1() # To compute the erosion height from all cells
sourceArea.chgele() # To substract the erosion height from original elevation value
sourceArea.diffusion() # To compute the diffusion for all the cells
if tt == 0:
if m2.pick == raise_slope_but.button and m2.release:
print "raise slope"
initial_slope += 0.1
if initial_slope > 0.6:
initial_slope = 0.6
print initial_slope
sourceArea.slope(dydz=initial_slope)
if m2.pick == lower_slope_but.button and m2.release:
print "lower slope"
initial_slope = initial_slope - 0.1
if initial_slope < 0.1:
initial_slope = 0.1
print initial_slope
sourceArea.slope(dydz=initial_slope)
[edit]
Gully erosion class (SAREA_6_tmp.py)
from visual import *
from random import *
from visual.graph import *
class area:
'''This class will create a source area'''
def __init__(self, nx, ny, dx, dy):
self.nx=nx # Number of rows
self.ny=ny # Number of columns
self.dx=dx # Unit size in x direction
self.dy=dy # Unit size in z direction
self.xmin = 0 # The area starting from zero x cord
self.ymax = nx*dx # The y cord is number of rows multiply by distance between two cells being top left corner as start of area
self.Elev=zeros((nx,ny), Float) # Initializing the array for elevation of all cells before providing slope
self.ElevT=zeros((nx*ny), Float) # Initializing the array for initial elevation of all cells
self.ElevS=zeros((nx,ny), Float) # Initializing the array for elevation of cells after providing slope
self.initElevS = zeros((nx,ny), Float) # Initializing the array as the initial elevation of cells after providing slope
self.land = [] # List for creating the different cells in the topography
seed(1000) # To make sure that every time same random numbers will be called by the program
for i in range(nx):
for j in range(ny):
self.land.append(box(length=dx, height=dy, # Appending the cell
width = (uniform(0.1,.5)),
color = (uniform(0.7,1),0.8,0.4)))
#up=(0,1.5,-0.26)))
self.land[-1].pos.z = self.land[-1].width / 2.0
self.land[-1].x=i*dx
self.land[-1].y=j*dy
self.Elev[i,j] = self.land[-1].width
#print 'the elevation', self.Elev
'''For better visual display changing the display characteristics'''
scene.center = (self.dx * (self.nx/2), self.dy * (self.ny/2), 0)
scene.forward = vector (1,0,0)
scene.up = vector (0,0,1)
scene.forward = vector(0,1,0)
scene.up = vector(0,0,1)
'''The following function is to provide slope to the flat topography'''
def slope(self, dxdz = 0, dydz = 0):
k=0
for i in self.land:
de = dxdz * i.pos.x
i.width = (uniform(0.1,.5)) + de
de = dydz * i.pos.y
i.width = i.width + de
i.pos.z = i.width / 2.0
self.ElevT[k] = i.width
k=k+1
hp=0
for i in range(self.nx):
for j in range(self.ny):
self.ElevS[i,j] = self.ElevT[hp]
self.initElevS [i,j]= self.ElevT[hp]
hp=hp+1
#print 'The elevation after slope'
#print self.ElevS
'''The following function is to give the value of number for given row and column'''
def ij_number(self,r=0,c=0):
Num = zeros((self.nx,self.ny))
num = 0
for i in range(self.nx):
for j in range(self.ny):
Num[i,j] = num
num=num+1
return Num[r,c]
'''The following function is to give the i and j value if number of that cell is provided'''
def number_ij(self,number):
T = number/self.nx
#print 'T', T
if T <=1:
i = 0
j = number-1
if (T > 1 and T <=2):
i = 1
j = number-self.nx-1
if (T >2 and T<=3):
i = 2
j = number - 2*self.nx-1
if (T >3 and T<=4):
i=3
j = number - 3*self.nx-1
if (T >4 and T<=5):
i = 4
j = number -4*self.nx-1
if (T>5 and T<=6):
i=5
j = number -5*self.nx-1
if (T>6 and T<=7):
i=6
j = number -6*self.nx-1
if (T>7 and T<=8):
i=7
j = number - 7*self.nx-1
if (T>8 and T<=9):
i=8
j = number -8*self.nx-1
if (T>9 and T<=10):
i=9
j = number - 9*self.nx-1
#print 'The i or Row is', i
#print 'The j or Column is', j
'''The following function is to provide the x and y co-ordinate for given R and C'''
def rc_xy(self,R,C):
self.x=self.xmin+self.dx/2.0+((self.dx)*C)
self.y=self.ymax-self.dx/2.0-((self.dx)*R)
return self.x, self.y
'''The following function is the subroutine for the function 'catchment' it calculates the
lowest value in that array of pixels '''
def LowRC(self):
Low = self.ElevS[self.PosR,self.PosC]
self.PosRL = self.PosR
self.PosCL = self.PosC
temp =0
for i in range(self.PosR-1,self.PosR+2):
for j in range(self.PosC-1,self.PosC+2):
if Low > self.ElevS[i,j]:
Low = self.ElevS[i,j]
self.PosRL = i
self.PosCL = j
return (self.PosRL, self.PosCL, Low)
'''The following function is also the subroutine for the function 'catchment' it sends the array for 3x3
surronding the given cell to function 'LowRC'''
def LowTrace(self,r,c): #
self.PosR = r
self.PosC = c
(LowRow, LowCol, Low) = self.LowRC()
self.cr_out[r,c] = self.ElevS[LowRow, LowCol]
self.c_out[r,c] = LowCol
self.r_out[r,c] = LowRow
'''The following function is to calculate the stream network i.e., the water flow through the system'''
def stream_network(self):
for i in range(1,self.nx-1):
for j in range(1,self.ny-1):
tr = curve(color=(uniform(0.1,1),0.1,0.7), radius = 0.13)
tmpr = i
tmpc = j
tr.append(pos=(i,j,self.ElevS[i,j]))
while not(self.c_out[tmpr,tmpc] == tmpc and self.r_out[tmpr,tmpc] == tmpr):
tmpc2 = self.c_out[tmpr, tmpc]
tmpr2 = self.r_out[tmpr, tmpc]
tmpc = tmpc2
tmpr = tmpr2
if not (tmpr == 0 and tmpc == 0):
tr.append(pos=(tmpr,tmpc,self.ElevS[tmpr,tmpc]))
'''The following function is to give the lowest elevation value adjacent to the given cell'''
def catchment(self,tt):
self.c_out = zeros((self.ny,self.nx))
self.r_out = zeros((self.nx,self.ny))
self.cr_out = zeros((self.nx,self.ny), Float)
for i in range(1,self.nx-1):
for j in range(1,self.ny-1):
self.LowTrace(i,j)
## print 'self.c_out', self.c_out
## print 'self.r_out', self.r_out
## print 'self.cr_out', self.cr_out
'''The following function calculates the watershed area i.e, number of the cells contributing to the given cell'''
def watershed(self,tt):
self.counter = ones((self.nx, self.ny))
seed(1000)
l_loop = 1
while l_loop == 1:
l_loop =0
for i in range(1,self.nx-1):
for j in range(1,self.ny-1):
temp = 1
for k in range(i-1,i+2):
for m in range(j-1,j+2):
if ((self.c_out[k,m] == j and self.r_out[k,m] == i) and not(k==i and m ==j)):
temp = temp +self.counter[k,m]
if temp <> self.counter[i,j]:
self.counter[i,j] = temp
l_loop = 1
(x,y) = self.rc_xy(i,j)
#print 'The catchment watershed area of area is',
#print self.counter
'''The following function calculates the local slope between the given cell and the lowest elevation cell'''
def localSlope(self):
self.Slope = zeros((self.nx,self.ny), Float)
for i in range(1,self.nx-1):
for j in range(1,self.ny-1):
self.Slope[i,j] = abs((self.cr_out[i,j]- self.ElevS[i,j])/self.dx)
#print 'the local slope is',
#print self.Slope
'''The following function calculates the sediment transport of all the cells'''
def sediTrans(self,tt):
#### '''The general sediment transport capacity is given as
#### Qs = (self.K* inv(self.ChiC)* self.ChiW* self.ChiT ^ self.P *
#### self.d ^ (1.5-self.P)* (1 - self.Tou* inv(self.ChiT) * (self.R*self.A)^ -self.MTou
#### * self.Slope^ -self.NTou) ^ self.P)* (self.R*self.A)^ -self.M * self.Slope^ -self.N'''
self.Qs = zeros((self.nx, self.ny), Float)
self.K = 20 # The calibration parameter
self.g = 9.81 # Acceleration due to gravity
self.P = 3 # Another calibration parameter
self.D = 0.003 # dominant grain size
self.TouCS = 0.045 # Dimensionless critical shear stress
self.R = .035 # Runoff rate in mm/hr
self.T = 0.25 # time for gully incision in hr
self.s = 2.65 # quartz specific gravity
self.RowW = 1000 # Water density
self.MTou = 0.6 # Exponent m
self.NTou = 0.8125 # Exponent n
self.M = 2.1 # Exponent M
self.N = 2.25 # Exponent N
Ks = 1.73
Kn = 0.08
C = 0.346
Ngc = 0.025
self.Qs = zeros((self.nx, self.ny), Float)
self.ChiC = pow (self.RowW,self.P) * pow((self.g* (self.s -1)), (self.P -0.5))
self.TouC = self.TouCS * (self.s-1) * self.D * self.RowW * self.g
self.ChiW = Ks*pow(Kn,0.375)*pow(C,-0.25)
self.ChiT = self.RowW*self.g*pow(C,0.75)*pow(Kn,-1.13)*pow(Ngc,1.5)
for i in range(1,self.nx-1):
for j in range(1,self.ny-1):
if ((self.counter[i,j] > 2) and (self.Slope[i,j] <> 0.0)):
self.A = self.counter[i,j]
self.Qs[i,j] = (self.K * (1/self.ChiC) * self.ChiW * pow(self.ChiT,self.P) * pow ( self.D, (1.5-self.P)) * pow((1 - self.TouC*(1/self.ChiT)* pow ((self.R * self.A), -self.MTou) * pow (self.Slope[i,j], -self.NTou)), self.P)) * pow ((self.R * self. A),self.M) * pow ((self.Slope[i,j]), self.N)
#self.land[self.ij_number(i,j)].color = color.red
#print 'The sediment capacity of different cells are',
#print self.Qs
'''The following function gives the erosion height of all cells depending upon the sediment capacity and time for gully incision'''
def eroht1(self):
AreaC = self.dx * self.dy
self.eroht = zeros((self.nx, self.ny), Float)
for i in range(self.nx):
for j in range(self.ny):
self.eroht[i,j] = (self.Qs[i,j]/AreaC)*self.T
#print 'erosion height of all cells', self.eroht
'''The following function computes the new elevation after substracting the erosion height'''
def chgele(self):
self.ChgElv = zeros((self.nx, self.ny), Float)
for i in range(1,self.nx-1):
for j in range(1,self.ny-1):
if self.eroht[i,j] > self.cr_out[i,j]:
self.eroht[i,j] = self.cr_out[i,j]
self.ChgElv[i,j] = self.ElevS[i,j]-self.eroht[i,j]
else:
self.ChgElv[i,j] = self.ElevS[i,j]-self.eroht[i,j]
#print 'change elevation', self.ChgElv
for i in range(1,self.nx-1):
for j in range(1,self.ny-1):
self.ElevS[i,j] = self.ChgElv[i,j]
#print 'final elev.',
#print self.ElevS
'''The following function is diffusion equation for all the cells'''
def diffusion(self):
'''Inflow is positive and outflow is negative'''
self.Qinx = zeros((self.nx, self.ny), Float)
self.Qoutx = zeros((self.nx, self.ny), Float)
self.Qiny = zeros((self.nx, self.ny), Float)
self.Qouty = zeros((self.nx, self.ny), Float)
self.deltaV = zeros((self.nx, self.ny), Float)
self.constant = .10
self.temp = zeros((self.nx, self.ny), Float)
for i in range(1, self.nx-1):
for j in range(1, self.ny-1):
self.Qinx[i,j] = self.constant* ((self.ElevS[i,j-1] - self.ElevS[i,j])/self.dx)
self.Qoutx[i,j] = -self.constant* ((self.ElevS[i,j] - self.ElevS[i,j+1])/self.dx)
self.Qiny[i,j] = self.constant* ((self.ElevS[i-1,j] - self.ElevS[i,j])/self.dx)
self.Qouty[i,j] = -self.constant* ((self.ElevS[i,j] - self.ElevS[i+1,j])/self.dx)
self.deltaV = self.Qinx + self.Qoutx + self.Qiny + self.Qouty
#print 'the deltaV', self.deltaV
for i in range(1, self.nx-1):
for j in range(1, self.ny-1):
self.temp[i,j] = self.deltaV[i,j]/ (self.dx*self.dy)
self.ElevS[i,j] = self.ElevS[i,j] + self.temp[i,j]
#print 'ele in diffusion'
#print self.ElevS
'''The following function is to update the elevation of topography'''
def update_topo(self):
## #print 1
for j in range(self.ny):
for i in range(self.nx):
self.land[self.ij_number(i,j)].width = self.ElevS[i,j]
self.land[self.ij_number(i,j)].pos.z = self.ElevS[i,j] / 2.0
'''The following function is to create new display for showing the final topography and erosion height '''
def newdisp(self,p):
if p ==1: # To show final eroded surface
tp = self.ElevS
else: # To show the diff. in elevation
tp = self.initElevS-self.ElevS
dispa = display ( x =0, y = 0, width= 640,height = 480)
for i in range(self.nx):
for j in range(self.ny):
self.land.append(box(length=self.dx, height=self.dy,
width = (uniform(0.1,.5))))
self.land[-1].x=i*self.dx
self.land[-1].y=j*self.dy
self.land[-1].width = tp[i,j]
self.land[-1].pos.z = self.land[-1].width / 2.0
if ((abs(self.temp[i,j]) > 0.09) and (p == 1)):
self.land[-1].color = color.blue
if (p ==0):
if tp[i,j] > 0:
self.land[-1].color = color.green
'''The following statements is for better display for the new display'''
dispa.center = (self.dx * (self.nx/2), self.dy * (self.ny/2), 0)
dispa.forward = vector (1,0,0)
dispa.up = vector (0,0,1)
dispa.forward = vector(0,1,0)
dispa.up = vector(0,0,1)
[edit]
Controls classes (uControls.py)
from visual import *
#see elipse.py for usage of slider
class slider:
def __init__(self, pos=vector(0,0,0),
min_val=0.0, max_val=1.0,
init_val=0.0,
stop_but= 0, stop_val=0.0,
scale=1.0,
title = "Control"):
'''
pos = position of the slider bar
min_val = minimum value on the slider bar
max_val = maximum value on the slider bar
init_val = initial value of slider bar (if beyond min_val or max_val is reset to min_val or max_val)
stop_but = set to 0 for no stop button and to 1 for a stop button (not yet in)
Note: The minimum y position of the slider is 0.0 and the maximum is self.scale
'''
self.pos = pos
self.min_val = min_val
self.max_val = max_val
self.stop_val = 0.0
self.scale = scale
self.value = init_val
self.title = title
self.bits = frame(pos=self.pos)
self.cyl = cylinder(frame=self.bits, radius=0.05*self.scale, axis=(0,self.scale,0), color=color.yellow)
self.but = sphere(frame=self.bits, pos=(0,self.get_y(), 0), radius = self.cyl.radius * 2.0, color = color.green)
self.title = label(frame=self.bits, pos=self.cyl.axis*1.2, text=self.title,height=10, box=0, xoffset=0, yoffset=0, space = 2)
self.val_label = label(frame=self.bits, pos=self.cyl.axis*-0.2, text=str(round(self.get_val(),1)),height=10, box=0, xoffset=0, yoffset=0, space = 2)
#adjust for boundary values
self.boundary_set()
def get_y(self):
return self.scale * (self.value - self.min_val) / (self.max_val - self.min_val)
def get_val(self):
return self.min_val + (self.but.pos.y / self.scale)*(self.max_val - self.min_val)
def boundary_set(self):
'''to keep values in bounds'''
if self.value < self.min_val:
self.value = self.min_val
self.but.pos.y = 0.0
if self.value > self.max_val:
self.value = self.max_val
self.but.pos.y = self.scale
def but_move(self, cursor_pos):
nx = 0.0 #self.bits.x
ny = cursor_pos.y
nz = 0.0 #self.bits.z
if ny < 0.0: #self.bits.y:
ny = 0.0 #self.bits.y
nv = self.min_val
elif ny > self.cyl.axis.y: #self.bits.y + self.cyl.axis.y:
ny = self.cyl.axis.y #self.bits.y + self.cyl.axis.y
nv = self.max_val
nv = self.min_val + ((ny) / self.scale)*(self.max_val - self.min_val)
print ny, nv
self.value = nv
self.val_label.text = str(round(self.get_val(),1))
return vector(nx, ny, nz)
class uSwitch:
def __init__(self, radius=1.0,init_val=0, pos=vector(0,0,0), title="Switch"):
self.title = title
self.value = init_val
self.bits = frame(pos=pos)
self.button = sphere(frame=self.bits, radius = radius, color = self.get_color())
self.title = label(frame=self.bits, pos=(0,radius*1.2,0),
text=self.title, height=10, box=0,
xoffset=0, yoffset=0, space = 2)
self.val_label = label(frame=self.bits, pos=(0,-radius*1.2,0),
text=self.get_val(),height=10, box=0, xoffset=0, yoffset=0, space = 2)
def switch(self):
if self.value == 1:
self.value = 0
#self.button.color = color.red
else:
self.value = 1
#self.button.color = color.blue
self.val_label.text = self.get_val()
self.button.color = self.get_color()
#print self.get_val()
def get_val(self):
if self.value == 1:
val = "On"
else:
val = "Off"
return val
def get_color(self):
if self.value == 1:
col = color.red
else:
col = color.blue
return col
