Code for gully erosion developed (modified version) (8th Nov. 2005)
From GeoMod
[edit]
Code for class 'area' created for gully erosion
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 teo cells being top left corner as start of area
self.Elev=zeros((nx,ny), Float)
self.ElevT=zeros((nx*ny), Float)
self.ElevS=zeros((nx,ny), Float)
self.land = []
seed(100)
for i in range(nx):
for j in range(ny):
self.land.append(box(length=dx, height=dy,
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
## colval = (self.land[-1].width-0.1)/(4-0.1)
## self.land[-1].color = (colval,colval,colval)
self.Elev[i,j] = self.land[-1].width
print 'the elevation', self.Elev
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)
def draax (self):
x_axis=cylinder(axis=(104,0,0), color=color.red, radius=2)
y_axis=cylinder(axis=(0,104,0), color=color.green, radius=2)
z_axis=cylinder(axis=(0,0,104), color=color.blue,radius=2)
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
p=0
for i in range(self.nx):
for j in range(self.ny):
self.ElevS[i,j] = self.ElevT[p]
p=p+1
print 'The elevation after slope'
print self.ElevS
## print "local Elevation = "
## print self.ElevS[1:4, 0:3]
return self.ElevS
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
## print 'the no of land pieces', Num
## print 'the number of r row and c col is', Num[r,c]
return Num[r,c]
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
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)
## print 'The x-cord for given Column, C is', C, self.x
## print 'The y-cord for given Row, R is', R, self.y
return self.x, self.y
def LowRC(self): # This function give array of pixels around the given pixel and lowest value in that array
Low = self.ElevS[self.PosR,self.PosC]
self.PosRL = self.PosR
self.PosCL = self.PosC
temp =0
## print 'The array of eight cell surronding the given cell'
## print self.ElevS[self.PosR-1:self.PosR+2,self.PosC-1:self.PosC+2]
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)
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
def catchment(self):
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
def watershed(self):
self.counter = ones((self.nx, self.ny))
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
## print self.counter
#rate(0.1)
print 'The catchment watershed area of area is',
print self.counter
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
def sediTrans(self):
#### '''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 = 1 # time for gully incision in hr
self.s = 2.65 # quartz specific gravity
self.RowW = 1000 # Specific gravity of water
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
## print 'here in sedi'
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)
## print "self.N", self.N
## print "self.ElevS[3,8]", self.Slope[3,8]
## print "pow ((self.ElevS[i,j]), self.N)", pow ((self.Slope[6,3]), self.N)
## print "pow ((self.R * self. A),self.M)", pow ((self.R * self. counter[0,0]+1),self.M)
## print "self.NTou", -self.NTou
## print "pow (self.ElevS[i,j], -self.NTou", pow (self.Slope[6,3], -self.NTou)
## print "pow ((self.R * self.A), -self.MTou)", pow ((self.R * self.counter[1,1]), -self.MTou)
## print "pow((1 - self.TouC*(1/self.ChiT)* pow ((self.R * self.A), -self.MTou) * pow (self.ElevS[i,j], -self.NTou)), self.P)) ",
## print pow((1 - self.TouC*(1/self.ChiT)* pow ((self.R * self.counter[0,0]+1), -self.MTou) * pow (self.Slope[6,3], -self.NTou)), self.P)
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)):
## print 'i,j', i, j
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
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
#return self.eroht
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.counter[i,j] > 3):
## self.ChgElv[i,j] = self.ElevS[i,j]+self.eroht[i,j]
## else:
if self.eroht[i,j] > self.ElevS[i,j]:
## print 'the eroht > elevs', self.eroht[i,j], self.ElevS[i,j]
self.eroht[i,j] = self.ElevS[i,j]
self.ChgElv[i,j] = self.ElevS[i,j]-self.eroht[i,j]
## print 'erosion ht ', self.eroht[6,3], self.eroht[i,j], self.ElevS[i,j]
## print 'chg elve after equating', self.ChgElv[i,j]
else:
## print 'eroht not greater'
## print self.eroht[i,j], self.ElevS[i,j]
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
def newdisp(self,t=1):
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))))
#color = (uniform(0.7,1),0.8,0.4)))
self.land[-1].x=i*self.dx
self.land[-1].y=j*self.dy
self.land[-1].width = self.ElevS[i,j]
self.land[-1].pos.z = self.land[-1].width / 2.0
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)
## self.land[self.ij_number(i,j)].color = color.green
[edit]
Code for main function calling this class is
from visual import *
from random import *
from SAREA_3 import *
from visual.graph import *
import math
# Define parameters
nx = 10
lx = 10.0
ny = 10
ly = 10.0
dx = lx/nx
dy = ly/ny
##print 'dx', dx
##print 'dy', dy
# creating land area as per the parameters defined
sourceArea = area (nx, ny, dx, dy)
# Translating the land to the center
#sourceArea.trans ( (nx-1)*lx/2, 0, (nz-1)*lz/2)
# drawing the axis
#sourceArea.draax()
sourceArea.slope(dydz=0.5)
# To get the land number from row and column
#sourceArea.ij_number(1,2)
# To get the value of i,j if the land number given
##sourceArea.number_ij(95.0) # Specify number as float quantity
#To get the contributing area of each pixel by noting the number of counter
tt = 15
while tt<>0:
sourceArea.catchment()
sourceArea.watershed()
sourceArea.localSlope()
sourceArea.sediTrans()
sourceArea.eroht1()
sourceArea.chgele()
if tt == 1:
sourceArea.newdisp(1)
tt= tt-1
--Ekta2 16:26, 6 Nov 2005 (CST)
