Final version of code of gully erosion
From GeoMod
from visual import *
from random import *
from SAREA_6_tmp import *
from visual.graph import *
import math
# 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
sourceArea.slope(dydz=0.5)
tt = 16 # Total number of times the code will execute
while tt<>0:
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 == 1):
sourceArea.newdisp(p=1) # To create the new display for the final topography
sourceArea.newdisp(p=0) # To create the new display for the total erosion from all cells
tt= tt-1
The code for class SAREA_6_tmp is as follows:
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 = i.width + 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):
for i in range(self.nx):
for j in range(self.ny):
self.land[self.ij_number(i,j)].width = self.ElevS[i,j]
self.land[self.ij_number(i,j)].pos.z = self.land[self.ij_number(i,j)].width / 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)
--Ekta2 15:34, 23 Nov 2005 (CST)
