Code for UM main campus - draft
From GeoMod
[edit]
This is the draft code for Nov.14th
The flood water looks too shallow, I'll try to figure out if it is true.
topography data:
water class:
from Numeric import *
from visual import *
from string import *
from random import uniform
from water_tile import *
'''Some thoughts and codes borrowed from Dr. Urbano's raster_map class. '''
class water:
def __init__(self, rasterfile, preciptation = 0.00000353, time = 20000): #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)-80 # standardization, relative elevation led the output looks better.
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))
self.elev = zeros((nrows, ncols,), Float)
self.Q = zeros((nrows, ncols,), Float)
self.P = preciptation
self.t = time
def slope1(self, i1, j1, i2, j2):
slope1 = (self.data[i1,j1] - self.data[i2,j2])/self.dx
return slope1
def slope2(self, i1, j1, i2, j2):
slope2 = (self.elev[i1,j1] - self.elev[i2,j2])/self.dx
return slope2
def d_tile_N(self,i,j):
if ((self.Q[i-1,j] > 0) and (sin(atan(self.slope1(i-1,j,i,j)))*cos(atan(self.slope1(i-1,j,i,j))) > 0)):
'''average horizontal vlocity,(0 + g*sin(A))/2), the total time,t1, to remove the water in [i-1,j]: dx/(0 + g*sin(A))/2)
the percentage of water moved: self.t/t1'''
## a = 9.81* sin(atan(self.slope1(i-1,j,i,j)))*cos(atan(self.slope1(i-1,j,i,j)))
## t1 = sqrt(2*dx/a)
## the ratio of water been moved: 1/t1
## #print v, t1
## d_tile_N = self.t/t1 * self.Q[i-1,j]
d_tile_N = 1/sqrt(2*self.dx/(9.81* sin(atan(self.slope1(i-1,j,i,j)))*cos(atan(self.slope1(i-1,j,i,j)))))* self.Q[i-1,j]
elif ((self.Q[i-1,j] > 0) and (sin(atan(self.slope1(i-1,j,i,j)))*cos(atan(self.slope1(i-1,j,i,j))) < 0)):
d_tile_N = -1/sqrt(-2*self.dx/(9.81* sin(atan(self.slope1(i-1,j,i,j)))*cos(atan(self.slope1(i-1,j,i,j)))))* self.Q[i-1,j]
else:
d_tile_N = 0
return d_tile_N
def d_tile_S(self,i,j):
if ((self.Q[i+1,j] > 0) and (sin(atan(self.slope1(i+1,j,i,j)))*cos(atan(self.slope1(i+1,j,i,j))) > 0)):
d_tile_S = 1/sqrt(2*self.dx/(9.81* sin(atan(self.slope1(i+1,j,i,j)))*cos(atan(self.slope1(i+1,j,i,j)))))* self.Q[i+1,j]
elif ((self.Q[i+1,j] > 0) and (sin(atan(self.slope1(i+1,j,i,j)))*cos(atan(self.slope1(i+1,j,i,j))) < 0)):
d_tile_S = -1/sqrt(-2*self.dx/(9.81* sin(atan(self.slope1(i+1,j,i,j)))*cos(atan(self.slope1(i+1,j,i,j)))))* self.Q[i+1,j]
else:
d_tile_S = 0
return d_tile_S
def d_tile_W(self,i,j):
## print("test"), sin(atan(self.slope1(i,j-1,i,j)), cos(atan(self.slope1(i,j-1,i,j))
if ((self.Q[i,j-1] > 0) and (sin(atan(self.slope1(i,j-1,i,j)))*cos(atan(self.slope1(i,j-1,i,j))) > 0)):
d_tile_W = 1/sqrt(2*self.dx/(9.81* sin(atan(self.slope1(i,j-1,i,j)))*cos(atan(self.slope1(i,j-1,i,j)))))* self.Q[i,j-1]
elif ((self.Q[i,j-1] > 0) and (sin(atan(self.slope1(i,j-1,i,j)))*cos(atan(self.slope1(i,j-1,i,j))) < 0)):
d_tile_W = -1/sqrt(-2*self.dx/(9.81* sin(atan(self.slope1(i,j-1,i,j)))*cos(atan(self.slope1(i,j-1,i,j)))))* self.Q[i,j-1]
else:
d_tile_W = 0
return d_tile_W
def d_tile_E(self,i,j):
if ((self.Q[i,j+1] > 0) and (sin(atan(self.slope1(i,j+1,i,j)))*cos(atan(self.slope1(i,j+1,i,j))) > 0)):
d_tile_E = 1/sqrt(2*self.dx/(9.81* sin(atan(self.slope1(i,j+1,i,j)))*cos(atan(self.slope1(i,j+1,i,j)))))* self.Q[i,j+1]
elif ((self.Q[i,j+1] > 0) and (sin(atan(self.slope1(i,j+1,i,j)))*cos(atan(self.slope1(i,j+1,i,j))) < 0)):
d_tile_E = -1/sqrt(-2*self.dx/(9.81* sin(atan(self.slope1(i,j+1,i,j)))*cos(atan(self.slope1(i,j+1,i,j)))))* self.Q[i,j+1]
else:
d_tile_E = 0
return d_tile_E
def d_Q(self,i,j):
#print self.d_tile_N(i,j), self.d_tile_S(i,j), self.d_tile_W(i,j), self.d_tile_E(i,j),self.P
d_Q = self.d_tile_N(i,j) + self.d_tile_S(i,j) + self.d_tile_W(i,j) + self.d_tile_E(i,j) \
+ self.P* pow(self.dx,2)
self.Q[i,j] = self.Q[i,j] + d_Q
d_Q = self.Q[i,j]
if d_Q > 0.0:
self.elev[i,j] = self.data[i,j] + d_Q / (2*pow(self.dx,2))
## else:
## d_Q = 0
## self.elev[i,j] = self.data[i,j]
return d_Q
def QTest(self,i,j,t):
time = 0
while time < self.t:
self.d_Q(i,j)
time = time + 1
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[imin:imax,jmin:jmax]
#draw in topography
tmax = max(max(tdata))
tmin = min(min(tdata))
flood = []
time = 0
while time < self.t:
for i in range(imin,imax):
for j in range(jmin,jmax):
#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,(self.d_Q(i,j) / pow(self.dx,2)* scale),self.dx,(uniform(0.5,1),0,0))
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[i,j]*scale
flood[-1].width = self.d_Q(i,j)
## 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[25,25]
time = time + 1
rate(100)
water_tile class:
from Numeric import *
from visual import *
from random import uniform
class water_tile(box):
def __init__(self,pos = (0,0,0),length = 1, width = 1, height = 1, color = (0,uniform(0,1),0)):
box.__init__(self, length=length, width=width , height=height, color=color)
self.elev = 0
def draw(self,r = 10, c =8):
for i in range(r):
for j in range(c):
print self.elev
topo class:
from Numeric import *
from string import *
from visual import *
'''Some thoughts and codes borrowed from Dr. Urbano's raster_map class. '''
class topo:
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)-80 # standardization, relative elevation led the output looks better.
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):
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 area(self):
x = self.dx * self.ncols
y = self.dx * self.nrows
area = x * y
return area
def rainfall(self, time = 60, r_rate = 0.0002117):
'''r_rate: mm/minute, time: minutes'''
Q = r_rate * time * self.area()
return Q
def low_points(self,arr):
r1 = 0
c1 = 0
r2 = 0
c2 = 1
inter = 0
low_point1 = arr[0,0]
low_point2 = arr[0,1]
if low_point1 > low_point2:
inter = low_point1
low_point1 = low_point2
low_point2 = inter
c1 = 1
c2 = 0
for i in range(0,3):
for j in range(0,3):
if not(i == 1 and j == 1):
if arr[i,j] < low_point2:
if arr[i,j] < low_point1:
low_point2 = low_point1
r2 = r1
c2 = c1
low_point1 = arr[i,j]
r1 = i
c1 = j
if arr[i,j] > low_point1:
low_point2 = arr[i,j]
r2 = i
c2 = j
print("results"),r1,c1,low_point1,r2,c2,low_point2
return r1,c1,low_point1,r2,c2,low_point2
def slope(self, i1, j1, i2, j2):
slope = (self.data[i1,j1] - self.data[i2,j2])/self.dx
return slope
def d_tile_N(self,i,j):
d_tile_N = sin(atan(self.slope(i-1,j,i,j)))* pow(self.dx,2) * 9.81
return d_tile_N
def d_tile_S(self,i,j):
d_tile_S = sin(atan(self.slope(i+1,j,i,j)))* pow(self.dx,2) * 9.81
return d_tile_S
def d_tile_W(self,i,j):
d_tile_W = sin(atan(self.slope(i,j-1,i,j)))* pow(self.dx,2) * 9.81
return d_tile_W
def d_tile_E(self,i,j):
d_tile_E = sin(atan(self.slope(i,j+1,i,j)))* pow(self.dx,2) * 9.81
return d_tile_E
main program:
from visual import *
from topo import *
from water import *
#import topography
UM_map = topo("umned3_2.txt")
UM_flood = water("umned3_2.txt")
#UM_map.draw_map(0, 139, 0,102, 10, 1)
UM_map.line_3d(0, 50, 0,50, 10)
print UM_map.area()
print UM_map.rainfall(50)
#local_matrix = UM_map.local_3x3(vector(56,30,10))
#UM_map.low_points(local_matrix)
##print UM_map.dx
##print UM_map.slope(2,2,2,3)
##print("atan2(-1):"),atan(-1)
##print("sin(atan2(-0.2345)):"),sin(atan(-1)), sin(-pi/4)
##print pow(2,3), pow(3,2)
UM_flood.draw_flood(20, 50, 20,50, 100)
##print("topo:"), UM_flood.data[45,46],UM_flood.data[44,46],UM_flood.data[46,46],UM_flood.data[45,45],UM_flood.data[45,47]
##print("slope:"), UM_flood.slope1(44,46,45,46),UM_flood.slope1(46,46,45,46),UM_flood.slope1(44,45,45,46),UM_flood.slope1(44,47,45,46)
##print("d_tile:"), UM_flood.d_tile_N(45,46),UM_flood.d_tile_S(45,46),UM_flood.d_tile_W(45,46),UM_flood.d_tile_E(45,46)
##print("d_water:"), UM_flood.d_water_N(45,46),UM_flood.d_water_S(45,46),UM_flood.d_water_W(45,46),UM_flood.d_water_E(45,46)
##print("Q:"), UM_flood.d_Q(45,46),UM_flood.d_Q(44,46),UM_flood.d_Q(47,46),UM_flood.d_Q(45,45),UM_flood.d_Q(45,47)
##print("uinit area:"), pow(UM_flood.dx,2)
##print("elev:"), UM_flood.elev[45,46],UM_flood.elev[44,46],UM_flood.elev[46,46],UM_flood.elev[45,45],UM_flood.elev[45,47]
##print ("d_tile_N(45,46)"),UM_flood.d_tile_N(45,46)
#UM_flood.QTest(45,46,20)
