10. Rising water level in one dimensional river - 14th Nov
From GeoMod
Here is my main program:
from visual import *
from raster_import import *
from Classes import *
#import topography
topog = raster_map("raster_ascii552.txt")
topog.data[0,1] = 500.0
topog.line_3d(0, 5, 0, 5, 0.01, 1)
#topog.draw_boundary()
topog.draw_map(0, 5, 0, 5, 0.01, 1)
#initialize water levels
topog.init_waterlevels(300)
while 1:
#add water to river
topog.wateradd(0.05)
#floodplain (scene.center, 1)
#topog.waterflow()
The program with functions and raster import goes:
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):
lowest_cell = 0
for i in range (self.ncols):
#print i, self.data[0,i]
if self.data[0,i] < self.data[0,lowest_cell]:
lowest_cell = i
#print "lowest_cell = ", lowest_cell
self.waterlevels[0,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 + lowest_cell* self.dx
jpos = self.ymax - 0 * self.dx
kpos = self.waterlevels[0,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,lowest_cell]*self.scale,
color = (0,0,1))
