Tracing a water droplet
From GeoMod
I thought that this would work... I retrieved some hints from Ekta's code but it didn't seem to do any good
Code for tracing a water droplet through the low points of a rasterfile
from Numeric import *
from string import *
from visual import *
from random import *
class raster_map:
def __init__(self, rasterfile): #datarr, xmin, ymax, dx, nrows, ncols):
'''rasterfile is the name given to the ArcGIS imported file'''
infile = open(rasterfile,"r")#r= read from that file if w then = write
ct = 0
line = infile.readline()
while line:
ct += 1
l = split(strip(line))# recognizes each group as a text group& puts into a list(l)
if ct == 1: #we are in the first line
ncols = int(l[1]) #puts it as an integer and l(1) is the second thing in the list
if ct == 2: #we are in the second line (count=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 == 6:
nodata = l[1]
if ct >= 7: #past line 7= topographic data
if ct == 7:
topo = zeros((nrows*ncols,), Float)#of rows and columns in size
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)
print 'The elevation at this position is',self.data[r,c]
return self.data[r,c]
def get_rc(self, pos):
'''pos = a vector value with x and y coordinates
eg: pos = vector(x_val, y_val)
to call use: .get_rc(vector(x_val,Y-val))'''
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]
self.imin=imin
self.imax=imax
self.jmin=jmin
self.jmax=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'''
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 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 raindrop (self,x_pos,y_pos,center=1):
x=x_pos
y=y_pos
z= self.get_val(vector(x,y))*self.scale
droplet= sphere (pos=(x,y,z), radius=999.0, color=color.blue)
if center == 1:
scene.center = (x_pos,y_pos,0)
def colNum (self, x_pos): #this will give the column number the drop is on
Dif = x_pos - self.xmin + self.dx/2
NumCol = int(Dif /self.dx)
self.NumCol=NumCol
print 'The column is', NumCol
self.PosC = self.NumCol
return NumCol
def rowNum (self, y_pos): # To give row number the drop is on
Dif = - y_pos + self.ymax + self.dx/2
NumRow = int(Dif /self.dx)
self.NumRow=NumRow
self.PosR = self.NumRow
self.NumCol=self.colNum(x_pos)
print 'The row is', NumRow
return NumRow
def rain_xy (self,R,C):
self.x=self.xmin-self.dx/2+((self.dx)*C)
self.y=self.ymax+self.dx/2-((self.dx)*R)
print 'The x-cord for given Column', C, self.x
print 'The y-cord for given Row', R, self.y
return (self.x, self.y)
def LowRC(self): # This function give array of pixels around the given pixel and lowest value in that array
Low = self.data[self.PosR,self.PosC]
print 'The array of eight cell surronding the given cell'
print self.data[self.PosR-1:self.PosR+2,self.PosC-1:self.PosC+2]
for i in range(self.PosR-1,self.PosR+2):
for j in range(self.PosC-1,self.PosC+2):
if Low > self.data[i,j]:
Low = self.data[i,j]
self.PosRL = i
self.PosCL = j
print 'The lowest value is', Low
return (self.PosRL, self.PosCL, Low)
def LowTrace(self): # Tracing where the drop is being led by the low points
tr = curve(color=color.red, radius = 200)
self.PosR = self.NumRow
self.PosC = self.NumCol
while ((self.PosC >= 0 and self.PosC <= self.ncols) and
(self.PosR >= 0 and self.PosR <= self.nrows)):
(LowRow, LowCol, Low) = self.LowRC()
if (LowRow == self.PosR and LowCol == self.PosC):
break
else:
(x, y) = self.xy(LowRow,LowCol)
tr.append(pos=(x,y,Low*self.scale))
self.PosR = LowRow
self.PosC = LowCol
print 'the low row', LowRow
print 'the low col', LowCol
scene.center = (x,y,Low)
f_topo = "rastert_asciito1.txt"
topo = raster_map(f_topo)
#start_x = 623721
#start_y = 3300295
#the array is called topo.data
#print topo.data
#topo.draw_boundary()
#topo.draw_map(100, 200, 100, 200, 50, 1)
#topo.line_3d(100,200,100,200,50,1)
topo.line_3d(scale=50)
start_x = 623721
start_y = 3300295
topo.rowNum(start_y)
topo.colNum(start_x)
'''Make the raindrop:'''
topo.raindrop(start_x, start_y)
'''Trace Where it goes:'''
topo.LowTrace()
#rc=topo.get_rc(vector(start_x, start_y))
#elev= topo.get_val(vector(start_x, start_y))
#print rc
#raindrop_r= rc[0]
#raindrop_c= rc[1]
#print 'row = ',raindrop_r, 'column =', raindrop_c
#print 'the start elevation is ',elev
