11. Final ID Model
From GeoMod
Finally my 1D model is working. I'm working on it to make it 2D...
Main Program
from visual import *
from raster_import3 import *
from Classes import *
#import topography
topog = raster_map("raster_ascii55.txt")
topog.data[0,1] = topog.data[1,0]= topog.data[1,4] = 8.0
topog.line_3d(0, 5, 0, 5, 0.01, 1)
#topog.draw_boundary()
topog.draw_map(0, 5, 0, 5, 0.1, 1)
#initialize water levels
topog.init_waterlevels(6)
#timestep
dt = 0.001
for i in range(200):
print "TIMESTEP = ", i
#add water to river (boundary condition)
topog.wateradd(0.1*dt)
#propogate water downstream (diffusion)
topog.diffuse2(dt)
And here goes the classes:
from Numeric import *
from string import *
from visual import *
class raster_map:
def __init__(self, rasterfile): #datarr, xmin, ymax, dx, nrows, ncols):
infile = open(rasterfile,"r")
ct = 0
line = infile.readline()
while line:
ct += 1
l = split(strip(line))
if ct == 1:
ncols = int(l[1])
if ct == 2:
nrows = int(l[1])
if ct == 3:
xmin = float(l[1])
if ct == 4:
ymin = float(l[1])
if ct == 5:
dx = float(l[1])
if ct == 7:
nodata = l[1]
if ct >= 7:
if ct == 7:
topo = zeros((nrows*ncols,), Float)
t_ct = -1
for i in l:
t_ct += 1
if i == nodata: #set no-data values to zero
topo[t_ct] = 0.0
else:
topo[t_ct] = float(i)
line = infile.readline()
infile.close
ymax = ymin + (dx * (nrows-1)) # top left corner
'''datarr is a 2d array '''
self.xmin = xmin
self.ymax = ymax
self.dx = dx
self.nrows = nrows
self.ncols = ncols
self.data = resize(topo, (nrows, ncols))
def get_val(self, pos):
r, c = self.get_rc(pos)
return self.data[r,c]
def get_rc(self, pos):
r = int(((self.ymax - pos.y + self.dx/2.0) ) /self.dx )
c = int(((pos.x - self.xmin + self.dx/2.0) ) /self.dx)
return (r, c)
def draw_map(self, imin=0, imax=-99999, jmin=0, jmax=-99999, scale = 1.0, center = 0):
''' draw a specific region of the raster map'''
self.scale = scale
default_val = -99999
if imax == default_val:
imax = self.nrows
if jmax == default_val:
jmax = self.ncols
#print "drawing map "
#print "imin, imax ", imin, imax
#print "jmin, jmax ", jmin, jmax
tdata = self.data[imin:imax,jmin:jmax]
#draw in topography
tmax = max(max(tdata))
tmin = min(min(tdata))
for i in range(imin,imax):
for j in range(jmin,jmax):
print i, j, self.data[i,j]
ipos = self.xmin + j* self.dx
jpos = self.ymax - i* self.dx
kpos = self.data[i,j] #self.get_val(vector(ipos,jpos,0))
colval = (kpos - tmin) / (tmax - tmin)
box(pos=(ipos,jpos, 0.5*kpos*scale),
length = self.dx, height = self.dx, width = self.data[i,j]*scale,
color = (colval,colval,colval))
#center scene on area
#print self.xmin + jmin* self.dx, self.ymax - imax* self.dx
if center == 1:
xc = self.xmin + (self.dx * ((jmin+jmax)/2))
yc = self.ymax - (self.dx * ((imin+imax)/2))
print "center", xc, yc
scene.center = vector(xc, yc, 0)
def draw_boundary(self):
bound = curve(color=color.blue)
bound.append(pos=(self.xmin, self.ymax,self.data[0,0]))
bound.append(pos=(self.xmin+self.dx*self.ncols, self.ymax,self.data[0,0]))
bound.append(pos=(self.xmin+self.dx*self.ncols,
self.ymax-self.dx*self.nrows,self.data[0,0]))
bound.append(pos=(self.xmin,
self.ymax-self.dx*self.nrows,self.data[0,0]))
bound.append(pos=(self.xmin, self.ymax,self.data[0,0]))
def line_3d(self, imin=0, imax=-99999, jmin=0, jmax=-99999, scale=1.0, center=1):
''' draw a specific region of the raster map as a wire mesh
imin - minimum column value to be drawn
imax - maximum column value
jmin - miinimum row
jmax - maximum row
scale - to add a vertical exageration
center - to center scene on area drawn if center == 1'''
self.scale = scale
default_val = -99999
if imax == default_val:
imax = self.nrows
if jmax == default_val:
jmax = self.ncols
print
tdata = self.data[imin:imax,jmin:jmax]
self.imin=imin
self.imax=imax
self.jmin=jmin
self.jmax=jmax
for i in range(imin, imax):
lin = curve(color=(0,1,0))
for j in range (jmin, jmax):
ipos = self.xmin + j* self.dx
jpos = self.ymax - i* self.dx
kpos = self.data[i,j]
lin.append(pos=(ipos,jpos,kpos*scale))
#print lin.pos[-1]
for j in range (jmin, jmax):
lin = curve(color=(0,1,0))
for i in range(imin, imax):
ipos = self.xmin + j* self.dx
jpos = self.ymax - i* self.dx
kpos = self.data[i,j]
lin.append(pos=(ipos,jpos,kpos*scale))
#center of scene of area
#print self.xmin + jmin* self.dx, self.ymax - imax* self.dx
if center == 1:
xc = self.xmin + (self.dx * ((jmin+jmax)/2))
yc = self.ymax - (self.dx * ((imin+imax)/2))
print "center", xc, yc
scene.center = vector(xc, yc, 0)
def local_3x3(self, pos):
#find the local 3x3 matrix of topography around the cell centered at pos
r, c = self.get_rc(pos)
#print "r,c=", r, c
return self.data[r-1:r+2, c-1:c+2]
def local_rc_3x3(self, r, c):
#find the local 3x3 matrix of topography around the cell centered at pos
#r, c = self.get_rc(pos)
#print "r,c=", r, c
return self.data[r-1:r+2, c-1:c+2]
def init_waterlevels(self, initial_elevation):
'''intitial_elevations is a single value: assuming initial water level is flat'''
self.waterlevels = float(initial_elevation) * ones([self.nrows, self.ncols])
self.waterdepths = self.waterlevels - self.data
print "waterlevels"
print self.waterlevels
print "waterdepths"
print self.waterdepths
# draw water
for i in range(self.nrows):
for j in range(self.ncols):
print i, j, self.waterlevels[i,j]
if self.waterlevels[i,j] > self.data[i,j]:
ipos = self.xmin + j* self.dx
jpos = self.ymax - i* self.dx
kpos = self.waterlevels[i,j] #self.get_val(vector(ipos,jpos,0))
box(pos=(ipos,jpos, 0.5*kpos*self.scale),
length = self.dx, height = self.dx, width = self.waterlevels[i,j]*self.scale,
color = (0,0,1))
def wateradd (self, quant):
self.lowest_cell = 0
for i in range (self.ncols):
#print i, self.data[0,i]
if self.data[0,i] < self.data[0,self.lowest_cell]:
self.lowest_cell = i
self.waterlevels[0,self.lowest_cell] += quant / (self.dx*self.dx)
print "quant", quant
print "dx", self.dx
#print "waterlevel", self.waterlevels[0,lowest_cell]
#draw change
ipos = self.xmin + self.lowest_cell* self.dx
jpos = self.ymax - 0 * self.dx
kpos = self.waterlevels[0,self.lowest_cell] #self.get_val(vector(ipos,jpos,0))
box(pos=(ipos,jpos, 0.5*kpos*self.scale),
length = self.dx, height = self.dx, width = self.waterlevels[0,self.lowest_cell]*self.scale,
color = (0,0,1))
def diffuse2(self, dt, mann=0.035):
j = 2
self.mann = mann
self.dt = dt
new_wl = self.waterlevels
for i in range(1, self.nrows-1):
#print " "
#print i
if self.waterlevels[i,j] > self.data[i,j]:
V_ym = 0.0
V_yp = 0.0
#For trailing face
if self.waterlevels[i,j]<self.waterlevels[i-1,j]:
slope = (self.waterlevels[i-1,j] - self.waterlevels[i,j]) / self.dx #calculating water surface elevation
self.ws_elev = self.waterdepths + self.data
hmax = self.hmax_define(float(self.ws_elev[i,j]), float(self.ws_elev[i-1,j]))
zmax = max(self.data[i,j], self.data[i-1,j])
#Call function equations
dh = self.equations(zmax, hmax, slope)
new_wl[i,j] = self.waterlevels[i,j] + dh
kpos = new_wl[i,j]
#box((pos = (i+(self.dx/2),j+(self.dx/2), 0.5*kpos*self.scale),
#length=self.dx, height=self.dx, width=new_wl[i,j]*self.scale,
#color=(0,0,1))
else:
slope = (self.waterlevels[i,j] - self.waterlevels[i-1,j]) / self.dx #calculating water surface elevation
self.ws_elev = self.waterdepths + self.data
self.ws_elev = self.waterdepths + self.data
hmax = self.hmax_define(float(self.ws_elev[i,j]), float(self.ws_elev[i-1,j]))
zmax = max(self.data[i,j], self.data[i-1,j])
#Call function equation
dh = self.equations(zmax, hmax, slope)
new_wl[i-1,j]=self.waterlevels[i-1,j] + dh
#For Forward Face
if self.waterlevels[i,j]<self.waterlevels[i+1,j]:
slope = (self.waterlevels[i+1,j] - self.waterlevels[i,j]) / self.dx #calculating water surface elevation
self.ws_elev = self.waterdepths + self.data
hmax = self.hmax_define(float(self.ws_elev[i,j]), float(self.ws_elev[i+1,j]))
zmax = max(self.data[i,j], self.data[i+1,j])
#call function equation
dh = self.equations(zmax, hmax, slope)
new_wl[i,j] = self.waterlevels[i,j] + dh
else:
slope = (self.waterlevels[i,j] - self.waterlevels[i+1,j]) / self.dx #calculating water surface elevation
self.ws_elev = self.waterdepths + self.data
hmax = self.hmax_define(float(self.ws_elev[i,j]), float(self.ws_elev[i+1,j]))
zmax = max(self.data[i,j], self.data[i+1,j])
#call function equation
dh = self.equations(zmax, hmax, slope)
new_wl[i+1,j]=self.waterlevels[i+1,j] + dh
print "New Waterlevel"
print new_wl
for i in range(self.nrows-1):
self.waterlevels[i,j] = new_wl[i,j]
self.waterlevels[i-1,j] = new_wl[i-1,j]
self.waterlevels[i+1,j] = new_wl[i+1,j]
def equations (self, zmax, hmax, slope):
#print "waterlevels"
#print self.waterlevels
#print "Data", self.data
diff = (hmax - zmax)
cs_area = diff * self.dx
w_perimeter = self.dx + (2*diff)
h_radius = (cs_area * self.dx ) / w_perimeter
flowrate_u = (cs_area * pow(h_radius, 2.0/3.0) * pow(abs(slope), 0.5)) / self.mann
V_ym = flowrate_u * self.dt #volume of water coming into cell
dh = V_ym / self.dx *self.dx
#print "slope = ", slope
#print "hmax = ", hmax
#print "zmax = ", zmax
#print "cs_area = ", cs_area
#print "w_perimeter = ", w_perimeter
#print "h_radius = ", h_radius
#print "flowrate_u = ", flowrate_u
#print "New Waterlevel"
#print new_wl
return dh
def hmax_define (self, a, b):
hmax = 0.0
if a > b:
hmax = a
else:
hmax = b
return hmax
