Code for Gravity Gradient Based Algorithm
From GeoMod
Basic logic of gravity gradient based algorithm:
Data for the project: Image:Umned3 2.txt
Code for class water11, gravity gradient based algorithm to modeling the flood inundation:
from Numeric import *
from visual import *
from string import *
from random import uniform
from visual.graph import * # import graphing features
from water_tile import *
'''Some thoughts and codes borrowed from Dr. Urbano's raster_map class. '''
'''This file is based on gravity gradient algorithm.'''
class water:
def __init__(self, rasterfile,imin=0, imax=-99999, jmin=0, jmax=-99999, preciptation = 0.00000353, time = 4300,t_step =0.2): #datarr, xmin, ymax, dx, nrows, ncols):
## '''a 'heavy rain' event is defined for this study as one that produces at least six inches in twelve hours or less
## at a single station, which is about the maximum amount of rain expected once every 100 years at a particular site
## in the Memphis CWA (Weather Bureau, 1961).
## 6 inches (152.399995 mm) /12 hr (12*3600 = 43,200 sec): 0.00353 mm/sc or 0.00000353 m/sc 0.2117 mm /minutes,
## 12.712 mm/hr '''
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
line = infile.readline()
infile.close
ymax = ymin + (dx * (nrows-1)) # top left corner
'''datarr is a 2d array '''
self.imin = imin
self.imax = imax
self.jmin = jmin
self.jmax = jmax
self.xmin = xmin
self.ymax = ymax
self.dx = dx
self.nrows = nrows
self.ncols = ncols
self.topo = resize(topo, (nrows, ncols))
self.elev = zeros((nrows, ncols,), Float)
self.o_elev = zeros((nrows, ncols,), Float)
for i in range(0,ncols-1):
for j in range(0,ncols-1):
self.elev[i,j] = self.topo[i,j]
self.o_elev[i,j] = self.elev[i,j]
self.Q = zeros((nrows, ncols,), Float)
self.o_Q = zeros((nrows, ncols,), Float)
self.depth = zeros((nrows, ncols,), Float)
self.P = preciptation
self.t = time
self.t_step = t_step
def get_lowest(self,imin,imax,jmin,jmax):
'''To get the lowest point in topo[imin:imax,jmin:jmax] and return the # of rows and columns '''
min_v = self.topo[imin,jmin]
r = imin
c = jmin
for i in range(imin,imax):
for j in range(jmin,jmax):
if self.topo[i,j] < min_v:
min_v = self.topo[i,j]
r = i
c = j
return r,c,min_v
def slope1(self, i1, j1, i2, j2):
'''get the slope of adjacent cells '''
slope1 = (self.topo[i1,j1] - self.topo[i2,j2])/self.dx
return slope1
def slope2(self, i1, j1, i2, j2):
'''get the slope of water in adjacent cells '''
slope2 = (self.o_elev[i1,j1] - self.o_elev[i2,j2])/self.dx
return slope2
## BASIC LOGIC OF GRAVITY GRADIENT ALGORITHM:
## 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_water_N = 1/t1 * self.Q[i-1,j]
## function d_water_N,d_water_S,d_water_W,d_water_E, and d_Q2 are defined to calculate water change
## based on gravity gradient
def d_water_N(self,i,j):
'''Determine the water communication with the adjacent cell on the North(i-1), positive: inflow; negative: outflow'''
d_water_N = 0
if j-1 >=self.jmin+1: #Do not calculate the boundary cell
if ((self.o_Q[i,j-1] > 0) and (sin(atan(self.slope2(i,j-1,i,j)))*cos(atan(self.slope2(i,j-1,i,j))) > 0)): # No water or no elevation difference,no flow
d_water_N = 1/sqrt(2*self.dx/(9.81* sin(atan(self.slope2(i,j-1,i,j)))*cos(atan(self.slope2(i,j-1,i,j)))))* self.o_Q[i,j-1] #inflow
if d_water_N > self.o_Q[i,j-1]: # the inflow can not exceed total amount availabe
d_water_N = self.o_Q[i,j-1]
elif ((self.o_Q[i,j] > 0) and (sin(atan(self.slope2(i,j-1,i,j)))*cos(atan(self.slope2(i,j-1,i,j))) < 0)): # elevation of water in adjacent cell lower, loose water
d_water_N = -1/sqrt(-2*self.dx/(9.81* sin(atan(self.slope2(i,j-1,i,j)))*cos(atan(self.slope2(i,j-1,i,j)))))* self.o_Q[i,j] #outflow
if d_water_N < -self.o_Q[i,j]: # the inflow can not exceed total amount availabe
d_water_N = -self.o_Q[i,j]
return d_water_N
def d_water_S(self,i,j):
'''Determine the water communication with the adjacent cell on the South(i+1), positive: inflow; negative: outflow'''
d_water_S = 0
if j+1 <= self.jmax-1: #Do not calculate the boundary cell
if ((self.o_Q[i,j+1] > 0) and (sin(atan(self.slope2(i,j+1,i,j)))*cos(atan(self.slope2(i,j+1,i,j))) > 0)): # No water or no elevation difference,no flow
d_water_S = 1/sqrt(2*self.dx/(9.81* sin(atan(self.slope2(i,j+1,i,j)))*cos(atan(self.slope2(i,j+1,i,j)))))* self.o_Q[i,j+1] #inflow
if d_water_S > self.o_Q[i,j+1]: # the inflow can not exceed total amount availabe
d_water_S = self.o_Q[i,j+1]
elif ((self.o_Q[i,j] > 0) and (sin(atan(self.slope2(i,j+1,i,j)))*cos(atan(self.slope2(i,j+1,i,j))) < 0)): # elevation of water in adjacent cell lower, loose water
d_water_S = -1/sqrt(-2*self.dx/(9.81* sin(atan(self.slope2(i,j+1,i,j)))*cos(atan(self.slope2(i,j+1,i,j)))))* self.o_Q[i,j] #outflow
if d_water_S < -self.o_Q[i,j]: # the inflow can not exceed total amount availabe
d_water_S = -self.o_Q[i,j]
return d_water_S
def d_water_W(self,i,j):
'''Determine the water communication with the adjacent cell on the West(j-1), positive: inflow; negative: outflow'''
d_water_W = 0
if i-1 >=self.imin+1: #Do not calculate the boundary cell
if ((self.o_Q[i-1,j] > 0) and (sin(atan(self.slope2(i-1,j,i,j)))*cos(atan(self.slope2(i-1,j,i,j))) > 0)): # No water or no elevation difference,no flow
d_water_W = 1/sqrt(2*self.dx/(9.81* sin(atan(self.slope2(i-1,j,i,j)))*cos(atan(self.slope2(i-1,j,i,j)))))* self.o_Q[i-1,j] #inflow
if d_water_W > self.o_Q[i-1,j]: # the inflow can not exceed total amount availabe
d_water_W = self.o_Q[i-1,j]
elif ((self.o_Q[i,j] > 0) and (sin(atan(self.slope2(i-1,j,i,j)))*cos(atan(self.slope2(i-1,j,i,j))) < 0)): # elevation of water in adjacent cell lower, loose water
d_water_W = -1/sqrt(-2*self.dx/(9.81* sin(atan(self.slope2(i-1,j,i,j)))*cos(atan(self.slope2(i-1,j,i,j)))))* self.o_Q[i,j] #outflow
if d_water_W < -self.o_Q[i,j]: # the inflow can not exceed total amount availabe
d_water_W = -self.o_Q[i,j]
return d_water_W
def d_water_E(self,i,j):
'''Determine the water communication with the adjacent cell on the East(j-1), positive: inflow; negative: outflow'''
d_water_E = 0
if i+1 <= self.imax-1: #Do not calculate the boundary cell
if ((self.o_Q[i+1,j] > 0) and (sin(atan(self.slope2(i+1,j,i,j)))*cos(atan(self.slope2(i+1,j,i,j))) > 0)): # No water or no elevation difference,no flow
d_water_E = 1/sqrt(2*self.dx/(9.81* sin(atan(self.slope2(i+1,j,i,j)))*cos(atan(self.slope2(i+1,j,i,j)))))* self.o_Q[i+1,j] #inflow
if d_water_E > self.o_Q[i+1,j]: # the inflow can not exceed total amount availabe
d_water_E = self.o_Q[i+1,j]
elif ((self.o_Q[i,j] > 0) and (sin(atan(self.slope2(i+1,j,i,j)))*cos(atan(self.slope2(i+1,j,i,j))) < 0)): # elevation of water in adjacent cell lower, loose water
d_water_E = -1/sqrt(-2*self.dx/(9.81* sin(atan(self.slope2(i+1,j,i,j)))*cos(atan(self.slope2(i+1,j,i,j)))))* self.o_Q[i,j] #outflow
if d_water_E < -self.o_Q[i,j]: #the outflow can not exceed total amount availabe
d_water_E = -self.o_Q[i,j]
return d_water_E
def d_Q(self,i,j):
'''calcuate the change of Q in every time step'''
d_Q = 0
d_Q += self.d_water_N(i,j) + self.d_water_S(i,j) + self.d_water_W(i,j) + self.d_water_E(i,j)
if d_Q < -self.Q[i,j]: # outflow can not exceed available amount
d_Q = -self.Q[i,j]
d_Q += self.P* pow(self.dx,2)*self.t_step #plus constant preciptation
return d_Q
def draw_flood(self,scale = 1.0, center = 0):
''' draw the flood water in the specific area of [imin:imax, jmin:jmax]'''
#create graph display
gdisp = gdisplay(x=0, y=0, width=480, height=350,title = "Depth vs. Time - Gravity Gradient Based",
ytitle="depth(m)(cyan:[r,c],orange:[r,c-1],red:[r,c-2])", xtitle="time(sc)")
# ymin=-50, ymax=110, xmax = 100)
Tcurve = gcurve(color=color.cyan) # a connected curve object
Tcurve2 = gcurve(color=color.orange)
Tcurve3 = gcurve(color=color.red)
self.scale = scale
default_val = -99999
if self.imax == default_val:
self.imax = self.nrows
if self.jmax == default_val:
self.jmax = self.ncols
## print "drawing map "
## print "imin, imax ", imin, imax
## print "jmin, jmax ", jmin, jmax
tdata = self.depth[self.imin:self.imax,self.jmin:self.jmax]
#draw in topography
tmax = max(max(tdata))
tmin = min(min(tdata))
colval = 0
flood = []
time = 0
print self.topo[self.imin:self.imax, self.jmin:self.jmax]
r,c,l = self.get_lowest(self.imin,self.imax,self.jmin,self.jmax)
print "Lowest point:"
print r,c,l
while time < self.t:
## if time > 1000:
## self.P=0
for i in range(self.imin,self.imax):
for j in range(self.jmin,self.jmax):
#idth = self.cal_depth(i,j)
#flood_water = water_tile((0,0,0),self.dx,self.dx,(self.d_Q(i,j) / pow(self.dx,2))*scale,(1,0,0))
dQ = self.d_Q(i,j)
if i == 135 and j==91:
print"Contribution of N,S,W,E:"
print self.d_water_N(i,j),self.d_water_S(i,j),self.d_water_W(i,j),self.d_water_E(i,j)
print("dQ of i,j,o_Q[i,j]:"),i,j, dQ, self.Q[i,j]
self.Q[i,j] = self.o_Q[i,j] +dQ
tp = self.topo[i,j]
## if i == 136 and j==92 :
## print ("i,j,elev[i,j]1,data[i,j]:"),i,j,self.elev[i,j],self.topo[i,j]
## self.elev[r,t] = tp + dQ/(2*pow(self.dx,2))
self.elev[i,j] = self.o_elev[i,j] + dQ/(2*pow(self.dx,2))
## if i == 136 and j==92 :
## print ("i,j,elev[r,t]2:"),i,j,self.elev[i,j]
if self.elev[i,j] < tp:
self.elev[i,j] = tp
## if i == 136 and j==92 :
## print ("i,j,elev[i]3:"),i,j,self.elev[i,j]
## print "total Q and total preceptation:"
## print "i,j:",i,j,self.Q[i,j],self.P*self.ncols*self.nrows*pow(self.dx,2)
self.depth[i,j] = self.elev[i,j] - self.topo[i,j]
#flood_water = water_tile((0,(self.elev[i,j]/2 + self.topo[i,j])*scale,0),self.dx,(self.elev[i,j]-self.topo[i,j])* scale,self.dx,(uniform(0.5,1),0,0))
flood_water = water_tile((0,0,0),self.dx,self.depth[i,j]* scale,self.dx,(colval,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]/2.0 + self.topo[i,j])*scale
flood[-1].z = self.elev[i,j]* scale
if tmax < 1:
colval = uniform(0.6,1)
else:
colval = (flood[-1].z - tmin) / (tmax - tmin)
print ("Maximum depth:")
print tmax,time
## if self.d_Q(i,j) > 0:
## flood[-1].width = self.d_Q(i,j) / pow(self.dx,2)* scale
## else:
## flood[-1].width = 0
flood[-1].depth = self.elev[i,j]-self.topo[i,j]
## if (i == 130 and j == 92):
#print("flood depth:"), flood[-1].depth
## flood[-1].width = (self.elev[i,j] - self.topo[i,j])*2*scale
## flood[-1].z = self.elev[i,j]*scale/2.0
## flood[-1].width = self.elev[i,j]*scale
for i in range(self.imin,self.imax):
for j in range(self.jmin,self.jmax):
self.o_elev[i,j] = self.elev[i,j]
self.o_Q[i,j] = self.Q[i,j]
if center == 1:
xc = self.xmin + (self.dx * ((self.jmin+self.jmax)/2))
yc = self.ymax - (self.dx * ((self.imin+self.imax)/2))
#print "center", xc, yc
scene.center = vector(xc, yc, 0)
print("elev[135,92],z, and topo:"), self.elev[135,91], self.elev[135,91]/2.0 + self.topo[135,91], self.topo[135,91]
print "Depth[135,92]:", self.depth[135,91]
time = time + self.t_step
print"tiem=:",time
Tcurve.plot(pos=(time, self.depth[r,c]))
Tcurve2.plot(pos=(time, self.depth[r,c-1]))
Tcurve3.plot(pos=(time, self.depth[r,c-2]))
## Tcurve.plot(pos=(time, self.depth[135,91]))
## Tcurve2.plot(pos=(time, self.depth[135,92]))
## Tcurve3.plot(pos=(time, self.depth[135,93]))
rate(100000000)
Code for class water_tile:
from Numeric import *
from visual import *
from random import uniform
'''A inheretant class from box, includes property of depth '''
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.depth = 0
Main program:
from visual import *
from topo import *
from water11 import *
#import topography
UM_map = topo("umned3_2.txt")
#UM_flood = water("umned3_2.txt",10,128,10,85)
UM_flood = water("umned3_2.txt",45,55,45,55,t_step =0.2)
#UM_map.draw_map(100, 138, 80,101, 10, 1)
UM_map.line_3d(0, 138, 0,101, 10, 1)
#UM_map.line_3d(0, 50, 0,50, 10)
#UM_map.draw_map(0, 138, 0,101, 10,1)
print UM_map.area()
print UM_map.rainfall(50)
##print "total preciptation:"
##print UM_flood.P *100*2000*pow(UM_flood.dx,2)
UM_flood.draw_flood(10, 1)
Code for class topo:
from Numeric import *
from string import *
from visual import *
'''A revised version of Dr. Urbano's raster_map class, add some functions. '''
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
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 get_lowest(self,imin,imax,jmin,jmax):
'''To get the lowest point in topo[imin:imax,jmin:jmax] and return the # of rows and columns '''
min_v = self.topo[imin,jmin]
r = imin
c = jmin
for i in range(imin,imax):
for j in range(jmin,jmax):
if self.topo[i,j] < min_v:
min_v = self.topo[i,j]
r = i
c = j
return r,c,min_v
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
