Test for 1 dimention
From GeoMod
The sample data for one dimention test for my algrithom: Umned3 4.txt
Main program:
from visual import *
#from topo import *
from water1 import *
UM_flood = water("umned3_4.txt")
UM_flood.QTest(14,200)
#UM_flood.draw_flood(0,1,0,14,10)
#UM_flood.cal()
class water1.py:
from Numeric import *
from visual import *
from string import *
from random import uniform
'''Some thoughts and codes borrowed from Dr. Urbano's raster_map class. '''
class water:
def __init__(self, rasterfile, preciptation = 0.00000353*1, time = 2000): #datarr, xmin, ymax, dx, nrows, ncols):
infile = open(rasterfile,"r")
ct = 0
nodata = -99999
line = infile.readline()
while line:
ct += 1
l = split(strip(line))
topo = zeros(15, 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)-86 # standardization
print i, topo[t_ct]
line = infile.readline()
infile.close
self.dx = 9.25
self.topo = topo
#self.elev = zeros(15, Float)
self.elev = self.topo
self.o_elev = self.elev
self.Q = zeros(15, Float)
self.o_Q = self.Q
self.P = preciptation
self.t = time
self.xmin = 0
self.ymax = 0
def slope1(self, j1, j2):
slope1 = (self.topo[j1] - self.topo[j2])/self.dx
return slope1
def slope2(self, j1, j2):
slope2 = (self.elev[j1] - self.elev[j2])/self.dx
return slope2
def right1(self,j):
right = 0
if j+1 <= len(self.topo):
## if j == 9 or j==1:
## print ("j,right, o_Q[j+1],sin,cos:"),j,(j+1), self.o_Q[j+1],sin(atan(self.slope1(j,j+1))),cos(atan(self.slope1(j,j+1)))
if ((self.o_Q[j+1] > 0) and (sin(atan(self.slope1(j+1,j)))*cos(atan(self.slope1(j+1,j))) > 0)):
right = 1/sqrt(2*self.dx/(9.81* sin(atan(self.slope1(j+1,j)))*cos(atan(self.slope1(j+1,j)))))* self.o_Q[j+1]
elif ((self.o_Q[j] > 0) and (sin(atan(self.slope1(j+1,j)))*cos(atan(self.slope1(j+1,j))) < 0)):
right = -1/sqrt(-2*self.dx/(9.81* sin(atan(self.slope1(j+1,j)))*cos(atan(self.slope1(j+1,j)))))* self.o_Q[j]
## if j == 9 or j==1:
## print("input of right1:"),right
## else:
## right = 0
return right
def left1(self,j):
left = 0
if (j-1)>= 0:
## if j == 9 or j==1:
## print ("j,left, o_Q[j-1],sin,cos:"),j,(j-1), self.o_Q[j-1], sin(atan(self.slope1(j-1,j))),cos(atan(self.slope1(j-1,j)))
if ((self.o_Q[j-1] > 0) and (sin(atan(self.slope1(j-1,j)))*cos(atan(self.slope1(j-1,j))) > 0)):
left = 1/sqrt(2*self.dx/(9.81* sin(atan(self.slope1(j-1,j)))*cos(atan(self.slope1(j-1,j)))))* self.o_Q[j-1]
elif ((self.o_Q[j] > 0) and (sin(atan(self.slope1(j-1,j)))*cos(atan(self.slope1(j-1,j))) < 0)):
left = -1/sqrt(-2*self.dx/(9.81* sin(atan(self.slope1(j-1,j)))*cos(atan(self.slope1(j-1,j)))))* self.o_Q[j]
## if j == 9 or j==1:
## print("input of left1:"), left
## else:
## left = 0
return left
def right2(self,j):
right = 0
if j+1 <= len(self.topo):
## if j == 9 or j==1:
## print ("j,right, o_Q[j+1],sin,cos:"),j,(j+1), self.o_Q[j+1],sin(atan(self.slope2(j,j+1))),cos(atan(self.slope2(j,j+1)))
if ((self.o_Q[j+1] > 0) and (sin(atan(self.slope2(j+1,j)))*cos(atan(self.slope2(j+1,j))) > 0)):
right = 1/sqrt(2*self.dx/(9.81* sin(atan(self.slope2(j+1,j)))*cos(atan(self.slope2(j+1,j)))))* self.o_Q[j+1]
elif ((self.o_Q[j] > 0) and (sin(atan(self.slope2(j+1,j)))*cos(atan(self.slope2(j+1,j))) < 0)):
right = -1/sqrt(-2*self.dx/(9.81* sin(atan(self.slope2(j+1,j)))*cos(atan(self.slope2(j+1,j)))))* self.o_Q[j]
## if j == 9 or j==1:
## print("input of right2:"),right
## else:
## right = 0
return right
def left2(self,j):
left = 0
if (j-1)>= 0:
## if j == 9 or j==1:
## print ("j,left, o_Q[j-1],sin,cos:"),j,(j-1), self.o_Q[j-1], sin(atan(self.slope2(j-1,j))),cos(atan(self.slope2(j-1,j)))
if ((self.o_Q[j-1] > 0) and (sin(atan(self.slope2(j-1,j)))*cos(atan(self.slope2(j-1,j))) > 0)):
left = 1/sqrt(2*self.dx/(9.81* sin(atan(self.slope2(j-1,j)))*cos(atan(self.slope2(j-1,j)))))* self.o_Q[j-1]
elif ((self.o_Q[j] > 0) and (sin(atan(self.slope2(j-1,j)))*cos(atan(self.slope2(j-1,j))) < 0)):
left = -1/sqrt(-2*self.dx/(9.81* sin(atan(self.slope2(j-1,j)))*cos(atan(self.slope2(j-1,j)))))* self.o_Q[j]
## if j == 9 or j==1:
## print("input of left2:"), left
## else:
## left = 0
return left
## def right(self,j):
## right = 0
## if j+1 <= len(self.data):
##
##
## print ("j,right, o_Q[j+1],sin,cos:"),j,(j+1), self.o_Q[j+1],sin(atan(self.slope2(j,j+1))),cos(atan(self.slope2(j,j+1)))
## if ((self.o_Q[j+1] > 0) and (sin(atan(self.slope2(j+1,j)))*cos(atan(self.slope2(j+1,j))) > 0)):
## right = 1/sqrt(2*self.dx/(9.81* sin(atan(self.slope2(j+1,j)))*cos(atan(self.slope2(j+1,j)))))* self.o_Q[j+1]
## elif ((self.o_Q[j] > 0) and (sin(atan(self.slope2(j+1,j)))*cos(atan(self.slope2(j+1,j))) < 0)):
## right = -1/sqrt(-2*self.dx/(9.81* sin(atan(self.slope2(j+1,j)))*cos(atan(self.slope2(j+1,j)))))* self.o_Q[j]
## print("input of right:"),right
## return right
## def left(self,j):
##
##
## left = 0
## if (j-1)>= 0:
## print ("j,left, o_Q[j-1],sin,cos:"),j,(j-1), self.o_Q[j-1], sin(atan(self.slope2(j-1,j))),cos(atan(self.slope2(j-1,j)))
## if ((self.o_Q[j-1] > 0) and (sin(atan(self.slope2(j-1,j)))*cos(atan(self.slope2(j-1,j))) > 0)):
## left = 1/sqrt(2*self.dx/(9.81* sin(atan(self.slope2(j-1,j)))*cos(atan(self.slope2(j-1,j)))))* self.o_Q[j-1]
## elif ((self.o_Q[j] > 0) and (sin(atan(self.slope2(j-1,j)))*cos(atan(self.slope2(j-1,j))) < 0)):
## left = -1/sqrt(-2*self.dx/(9.81* sin(atan(self.slope2(j-1,j)))*cos(atan(self.slope2(j-1,j)))))* self.o_Q[j]
## print("input of left:"), left
## return left
def d_Q(self,j):
#print ("left, right,p:"),self.left(j),self.right(j),self.P * pow(self.dx,2)
#d_Q = self.left1(j) + self.right1(j) + self.P * pow(self.dx,2)
d_Q = self.left2(j) + self.right2(j) + self.P * pow(self.dx,2)
#d_Q = self.left1(j) + self.right1(j) + self.left2(j) + self.right2(j)+ self.P * pow(self.dx,2)
#print("d_Q"), d_Q
## if d_Q >0.0:
##
## self.elev[j] += d_Q/(2*pow(self.dx,2))
return d_Q
def cal_depth(self,j):
dQ = self.d_Q(j)
self.Q[j] += dQ
self.o_Q[j] = self.Q[j]
if self.Q[j] > 0:
cal_depth = self.Q[j]/(2*pow(self.dx,2))
else:
cal_depth = 0
self.elev[j] = self.data[j] + cal_depth
## if j==11:
## print self.d_Q(j)
return cal_depth
def QTest(self,j,t):
#sphere()
#rate(100)
time = 0
table = []
while time < t:
#rate (100)
#print ("time ="),time
for i in range(0,j):
dQ = self.d_Q(i)
## if i == 9 or j==1:
## print("dQ of j:"),i, dQ
self.Q[i] += dQ
tp = self.topo[i]
if i == 8 or i==9 :
print ("i,elev[i]1,data[i]:"),i,self.elev[i],self.topo[i]
## self.elev[i] = tp + dQ/(2*pow(self.dx,2))
self.elev[i] = self.o_elev[i] + dQ/(2*pow(self.dx,2))
if i == 8 or i==9 :
print ("i,elev[i]2:"),i,self.elev[i]
if self.elev[i] < tp:
self.elev[i] = tp
if i == 8 or i==9 :
print ("i,elev[i]3:"),i,self.elev[i]
self.o_elev[i] = self.elev[i]
table.append(box(pos=(i-6,self.elev[i]/2,0),longth =2, width = 1, height =(self.elev[i]-self.topo[i])*100, color = (uniform(0.5,1),0,0)))
## if i == 8 or i==9 :
## print("i, Q[i], data[i], elev[i], depth:"), i, self.Q[i], self.topo[i], self.elev[i], self.elev[i] - self.topo[i]
## for i in range(0,j):
## print("Q of i:"),i,":", self.Q[i], self.elev[i], self.topo[i], self.elev[i]-self.topo[i]
time += 1
def cal(self):
tiles = []
for t in range(self.t):
for i in range(1,len(self.data)-1):
print i
self.elev[i] = (pow(self.o_elev[i-1],2) - 2* pow(self.o_elev[i],2) + pow(self.o_elev[i+1],2)) \
* (1/2*pow(self.dx,2)) + self.o_elev[i]
print self.elev[i]
tiles.append(box(pos =(i*self.dx,self.elev[i]/2.0,0), length = self.dx, height = self.elev[i], \
color = (uniform(0,1),0,0)))
tiles[-1].pos.y = self.elev[i]/2.0
tiles[-1].height = self.elev[i]
self.o_elev = self.elev
rate(10)
def draw_flood(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
#draw in topography
## tmax = max(max(tdata))
## tmin = min(min(tdata))
flood = []
time = 0
while time < self.t:
rate(50)
for i in range(imin,imax):
for j in range(jmin,jmax):
#print "j", j
width = self.cal_depth(j)
#print width
#flood_water = water_tile((0,0,0),self.dx,self.dx,(self.d_Q(i,j) / pow(self.dx,2))*scale,(1,0,0))
flood_water = water_tile((0,0,0),self.dx,width* scale,self.dx,(uniform(0,1),1,1))
flood.append(flood_water)
flood[-1].x = self.xmin + j* self.dx
flood[-1].y = self.ymax - i* self.dx
flood[-1].z = self.elev[j]*scale
## if self.Q[j] > 0:
## flood[-1].width = self.Q[j] / pow(self.dx,2)* scale
## else:
## flood[-1].width = 0
flood[-1].depth = width
## if j==13:
## print self.Q[j]
#print("j, flood depth:"),j, flood[-1].depth
## flood[-1].width = (self.elev[i,j] - self.data[i,j])*2*scale
## flood[-1].z = self.elev[i,j]*scale/2.0
## flood[-1].width = self.elev[i,j]*scale
#print self.elev[i,j]
#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)
#print self.elev[130,92]
time = time + 1
#print time
