The MAIN PROGRAM
From GeoMod
It's a bit messy but I hope you can follow!
from Numeric import *
from string import *
from visual import *
from raster_mapd import *
#DRAWING THE CHOSEN DEM TOPO MAP#
topo = raster_import_ryan("file.txt", dx=10.0)
topo.line_3d(scale=30,center=1)
#ball1=sphere()
#DEFINING VARIABLES#
starting_pos = vector(0,0,0)
starting_col= 0
starting_row= 0
highest_pos= vector(0,0,0)
lowest_pos= vector(0,0,0)
basin_arr=[]
nrows= topo.nrows
ncols= topo.ncols
#print nrows
#print topo.nrows, '=number of rows'
starting_pos.z=0 #default
starting_pos.x =topo.get_x(starting_row)
starting_pos.y = topo.get_y(starting_col)
starting_pos.z= topo.get_val(starting_pos)
#print 'Starting position', starting_pos
#print data
'''A= zeros((nrows,ncols))
for i in range (0,nrows) : #filling in array from top down
#print i
for j in range (0,ncols):
x=A[i,j]
#finding the highest # in the Array A
if x > high_num :
high_num=x
print 'the high number' ,x, 'came from row' ,i+1, 'and column' ,j+1
print A
print high_num'''
#GETTING THE HIGHEST POSITION IN THE CHOSEN TOPOGRAPHY#
(highest_r, highest_c, highest_elev) = topo.get_highest(starting_pos)
#print highest_r , '=highest row', highest_c, '= highest column'
#print 'the highest elevation is', highest_elev
highest_pos.z=highest_elev
highest_pos.x =topo.get_x(highest_c)
highest_pos.y = topo.get_y(highest_r)
#print highest_pos, '=highest position'
#print "-"
highpoint_topo= topo.local_3x3(highest_pos)
#print highpoint_topo
scale=30
#rod = cylinder(pos=highest_pos, axis=(5*scale,0,0), radius=1*scale)
topo.draw_post(highest_pos)
new_pos= highest_pos
new_c = highest_c
new_r = highest_r
#DEFINING THE DRAINAGE BASIN
#top_curve= curve(color=color.red, pos=new_pos, radius=5)
topo.draw_curve(new_pos)
#ball=sphere()
#ball=[]
#VARIABLES
default_elev = 0.0
local_high=0.0
old_local_high= 0.0
old_pos = new_pos
left_basin=[]
left_basin.append(highest_pos)
ct_left = 1
t=0
w=0
new_i = 1
new_j=1
local_r=1
local_c=1
#while (t<> 40):
# t=t+1
oldr = []
oldr.append(new_r)
oldc = []
oldc.append(new_c)
#while (new_pos.z <= old_pos.z):
while (local_r>0 and local_r< (nrows-1))and (local_c>0 and local_c< (ncols-1)):
#t=t+1
w=w+1
#highpoint_topo= topo.local_3x3(new_pos)
#local_arr = topo.data[oldr[-1]-1:oldr[-1]+2, oldc[-1]-1:oldc[-1]+2]
#print local_arr
local_high= default_elev
for i in range (oldr[-1]-1,oldr[-1]+2):
for j in range (oldc[-1]-1,oldc[-1]):
#print i, j, topo.data[i,j]
if not (i==oldr[-1] and j==oldr[-1]):
if (topo.data[i,j] > local_high):
dont_do = 0
#print ""
for n, t in enumerate(oldr):
#print n, oldr[n], oldc[n]
if (oldr[n] == i and oldc[n] == j): #dont do it
dont_do = 1
if dont_do == 0:
local_r= i
local_c= j
local_high= topo.data[i,j]
oldr.append(local_r)
oldc.append(local_c)
new_pos= vector( topo.get_x(local_c),topo.get_y(local_r), local_high)
left_basin.append(new_pos)
ct_left=ct_left+1
#ball.append(sphere(pos=vector(0,0,0),radius=6))
#print 'new position', new_pos, t
t=t+1
topo.append_curve(new_pos)
#print left_basin, 'left_basin positions'
#print 'ct for left basin', ct_left
#RIGHT SIDE OF DRAINAGE BASIN
topo.draw_curve(highest_pos)
new_c2 = highest_c
new_r2 = highest_r
#VARIABLES
default_elev = 0.0
local_high=0.0
old_local_high2= 0.0
new_pos=highest_pos
old_pos = new_pos
right_basin=[]
right_basin.append(highest_pos)
ct_right = 1
ct = 0
t2=0
w2=0
new_i2 = 1
new_j2=1
local_r=1
local_c=1
oldr2 = []
oldr2.append(new_r2)
oldc2 = []
oldc2.append(new_c2)
while (local_r>0 and local_r< (nrows-1))and (local_c>0 and local_c< (ncols-1)):
t2=t2+1
w2=w2+1
#highpoint_topo= topo.local_3x3(new_pos)
#local_arr = topo.data[oldr[-1]-1:oldr[-1]+2, oldc[-1]-1:oldc[-1]+2]
#print local_arr
local_high= default_elev
for i in range (oldr2[-1]+1,oldr2[-1]+2):
for j in range (oldc2[-1],oldc2[-1]+2):
#print i, j, topo.data[i,j]
if not (i==oldr2[-1] and j==oldr2[-1]):
if((topo.data[i,j])>(local_high)):
dont_do = 0
#print ""
for n, t2 in enumerate(oldr2):
#print n, oldr[n], oldc[n]
if (oldr2[n] == i and oldc2[n] == j): #dont do it
dont_do = 1
if dont_do == 0:
local_r= i
local_c= j
local_high= topo.data[i,j]
oldr2.append(local_r)
oldc2.append(local_c)
new_pos= vector( topo.get_x(local_c),topo.get_y(local_r), local_high)
right_basin.append(new_pos)
ct_right=ct_right+1
#print 'new position', new_pos, t2
topo.append_curve(new_pos)
#print right_basin, 'right_basin positions'
#print 'ct for left basin', ct_left, 'ct for right basin', ct_right
#FINDING AND DRAWING THE MIDPOINT LINE
topo.draw_curve(highest_pos)
mid_dis=[]
midpoint_pos=[]
ct_mid=-1
for i in range (1,ct_left):
dis_mid= topo.mid_dis(left_basin[i], right_basin[i])
mid_dis.append (dis_mid)
pos_midpoint= topo.midpoint (left_basin[i], right_basin[i])
midpoint_pos.append(pos_midpoint)
ct_mid= ct_mid+1
#drawing the horizontal and post line between points
topo.append_curve(midpoint_pos[i-1])
#topo.draw_post(midpoint_pos[i-1])
##########TESTING FOR LOWER ROW BOUNDARY###############
#print"-"
#print midpoint_pos[ct_left-2],'position of the bottom of midpoint line'
(low_row,low_col)= topo.get_rc(midpoint_pos[ct_left-2])
#print 'row', low_row, 'column',low_col
test_elev= topo.get_val(midpoint_pos[ct_left-2])
low_start=vector(midpoint_pos[ct_left-2].x,midpoint_pos[ct_left-2].y,test_elev)
topo.draw_post(low_start)
#print"-"
#######################################################
#THE BOTTOM BOUDARY
topo.draw_curve(left_basin[ct_left-2])
topo.append_curve(right_basin[ct_left-2])
#FINDING THE lOWEST POINTS WHICH ARE ASSUMED TO BE THE RIVER CHANNEL
default_elev2=9999
local_low=9999
topo.draw_curve(low_start)
new_c3 = low_col
new_r3 = low_row
new_pos=low_start
river_test=[]
river_test.append(low_start)
ct_river=0
t3=0
local_r=highest_r+2
local_c=1
oldr3=[]
oldr3.append(new_r3)
oldc3=[]
oldc3.append(new_c3)
##### FINDING THE RIVER ########################
while (local_r>=highest_r+2 and local_r< (nrows-1))and (local_c>0 and local_c< (ncols-1)):
#while (t3<=20):
t3=t3+1
#w3=w3+1
#highpoint_topo= topo.local_3x3(new_pos)
#local_arr = topo.data[oldr[-1]-1:oldr[-1]+2, oldc[-1]-1:oldc[-1]+2]
#print local_arr
local_low= default_elev2
for i in range (oldr3[-1]-1,oldr3[-1]):
for j in range (oldc3[-1]-1,oldc3[-1]+2):
#print i, j, topo.data[i,j]
if not (i==oldr3[-1] and j==oldr3[-1]):
if (topo.data[i,j] < local_low):
dont_do = 0
#print ""
for n, t3 in enumerate(oldr3):
#print n, oldr3[n], oldc3[n]
if (oldr3[n] == i and oldc3[n] == j): #dont do it
dont_do = 1
if dont_do == 0:
local_r= i
local_c= j
local_low= topo.data[i,j]
oldr3.append(local_r)
oldc3.append(local_c)
new_pos= vector( topo.get_x(local_c),topo.get_y(local_r), local_low)
river_test.append(new_pos)
ct_river=ct_river+1
#ball.append(sphere(pos=vector(0,0,0),radius=6))
#print 'new position', new_pos, t
topo.append_curve(new_pos)
#print river_test, 'initial river positions'
#print 'ct for the river array', ct_river
river=[]
#river.append(river_test[ct_river])
#print river
y=ct_river
while (y<=ct_river) and (y>=0):
river.append(river_test[y])
y=y-1
#print river
'''midpoint_pos[] and river[] and the left_basin[]... are the arrays being used'''
''' ct_mid, ct_river, ct_left'''
#print ct_mid,'ct_mid'
#print ct_river,'ct_river'
#print ct_left,'ct_left'
Da_array=[]
dir_vec=[]
i=0
j=15
ct_da=-1
ct_dir=-1
topo.draw_post(midpoint_pos[12])
topo.draw_post(left_basin[18])
########### FINDING DA (DISTANCE FROM RIVER TO MIDLINE) ###################
while ((i>=0) and (i<=ct_river)) and ((j>=15) and (j<=ct_river+15)):
Da= topo.distance(river[i], midpoint_pos[j])
Da_array.append(Da)
dir_=norm(river[i]-midpoint_pos[j])
dir_vec.append(dir_)
i=i+2
j=j+2
ct_da=ct_da+1
ct_dir=ct_dir+1
#print Da_array
#print ct_da, 'ct_da'
########### FINDING DD (DISTANCE FROM LEFT BOUNDARY TO MIDLINE) ############
k=15
l=18
Dd_array=[]
mag_ratio=[]
ct_dd=-1
while ((k>=15) and (k<=ct_river+15)) and ((l>=18) and (l<=ct_river+18)):
Dd= topo.distance(left_basin[l], midpoint_pos[k])
Dd_array.append(Dd)
k=k+2
l=l+2
ct_dd=ct_dd+1
#print Dd_array
#print ct_dd, 'ct_dd'
#print"-"
#print dir_vec
############### MAGNITUDE RATIO Da/Dd ######################
print "-"
for i in range (0,ct_da):
data_ratio= Da_array[i]/Dd_array[i]
mag_ratio.append(data_ratio)
print mag_ratio[i]
print '-'
#print mag_ratio
#print'-'
#print new_pos
############### Create a new screen window #############
scene2 = display( title = "Act IV, Scene 2" )
for i in range (0,ct_da):
next_pos=mag_ratio[i]*dir_vec[i]
vector_curve =curve( display = scene2, pos=[(0,0,0),(next_pos)], color=color.green)
#vector_curve.append (display = scene2,pos= next_pos, color=color.blue)
