Chertman.py
From GeoMod
To run the code, once you have VPython installed, you can either.
- Download this file: chert_man_24-rm.py
- or, copy and paste the following code to a file (named "chert_man_24-rm.py" for example) using your favorite text editor.
Other necessary files:
- Media:Rastert asciito1.txt
- Media:rechecked_xy_transf.txt
- Media:streams_transformed2.txt
- raster_map.py
'''
Chertman the model
Designed by
Dr. Lensyl Urbano
Dr. Stephen Cole
Department of Earth Sciences
University of Memphis
Published on web 2005
Free for non-commercial use
Help and support: Look us up at the University of Memphis website
'''
from visual import *
from random import *
from string import *
from visual.graph import * # import graphing features
from datetime import *
from Numeric import *
from Scientific.Statistics import *
#from Scientific.Geometry import Vector
from sys import *
import os
from os.path import basename
from raster_map import *
class outcrop:
def __init__(self, id, x, y, l_draw = 0, crange = 100.0, col = vector(1,0,1), z = 0.0,
Azero = 0.0, halflife = 9999.0):
self.x = float(x)
self.y = float(y)
self.id = id
self.Azero = Azero #Depth of exponential attractor function (A0)(see exponential_attractors.xls)
self.halflife = halflife #Half life of exponential attractor function (tau) (see exponential_attractors.xls)
if l_draw == 1:
#print type(self.x), type(self.y), type(z), type(crange), type(col)
self.ball = sphere(pos=(self.x, self.y, z), radius=crange, color=col)
def elevate(self, raster_map, z_scale = 1.0): #set the x value for the outcrop position based on a dem
self.ball.pos.z = raster_map.get_val(self.ball.pos) * z_scale
def elevate(pos, raster_map, z_scale = 1.0): #set the x value for the outcrop position based on a dem
pos.z = raster_map.get_val(pos) * z_scale
return pos
class stream:
def __init__(self, name):
self.x = [] #list of x locations along line
self.y = [] #list of y locations along line
self.z = [] #list of y locations along line
self.size = [] #stream size
self.name = name #name of stream
def cross(self): #this method will determine which stream reaches are crossable or not
self.crossable = equal(self.size, 1) # a value of 1 means NOT CROSSABLE
def elevate(self, raster_map, z_scale = 1.0):
#print "self.x = ", self.x
for i in range(len(self.x)):
self.z[i] = raster_map.get_val(vector(self.x[i], self.y[i], 0.0)) * z_scale
## print "stream = ", self.name
## print "crossable = ", self.crossable
class hit_counter:
def __init__(self, id):
self.id = id
self.ct = 0 #number of times outcrop type has been encountered
self.mass = 0.0 #mass of outcrop-type material in bag
self.avmass = 0.0 #average mass of material in bag after walk
self.stmass = [] #record of mass at the end of every walk
self.stdev = 0.0 #standard deviation of masses at end of walks
self.min = 0.0 #minimum of masses at end of walks
self.max = 0.0 #maximum of masses at end of walks
self.median = 0.0 #median of masses at end of walks
self.skew = 0.0 #skewness of masses at end of walks
self.atr_rate = 0.0 #attrition rate of type of material
self.prob_pickup = 1.0 #probability of picking up this type of material
class attractor:
def __init__(self, id, Azero = 0.0, halflife = 99999.0, isactive = 1.0):
self.id = id
self.Azero = Azero
self.halflife = halflife
self.isactive = isactive
class cherttype:
def __init__(self, id, Azero = 0.0, halflife = 9999999.0, attractive = 0.0, atr_rate = 0.0,
prob_pickup = 1.0):
self.id = id
self.Azero = Azero
self.halflife = halflife
self.attractive = attractive #if attractors are active is active
self.atr_rate = atr_rate
self.prob_pickup = prob_pickup
self.x = [] #array giving the x, y and z coordinates of each outcrop in type
self.y = []
self.z = []
def makearrays(self):
self.x = array(self.x)
self.y = array(self.y)
self.z = array(self.z)
def take_step(nhits,ndirs,step_size,oray,hit_ct,outcrops,load_add, max_load, total_mass):
#Take step
#print "newpos = ", newpos
#print newpos[0], newpos[1], newpos[2]
dx = newpos[0] - oray[:,1]
dy = newpos[1] - oray[:,2]
rad = sqrt(dx * dx + dy * dy) #calculates the distance to each point (outcrop)
##print rad
rad_intercept = less_equal(rad, rng) #determine if the distance to any point is less than the range of influence (rng)
rad_list = rad_intercept.tolist()
#print rad_intercept
hitlist = sort(nonzero(rad_intercept.flat)).tolist() # i,j encoded as i*Natoms+j
#print "hitlist, len(hitlist) = ", hitlist, len(hitlist)
nhits = nhits + len(hitlist)
# determine the id of outcrops encountered
if len(hitlist) > 0:
#print 'initial hitlist = ', hitlist
# calculate if material is picked up given the probability of pickup (prob_pickup)
newhitlist = []
for j in hitlist:
for k in hit_ct:
if outcrops[j].id == k.id:
if random() <= k.prob_pickup: #then pickup
newhitlist.append(j)
hitlist = newhitlist
#print 'new hitlist = ', hitlist
if len(hitlist) > 0:
# If multiple sites are selected with each step then add the same amount for each site picked
# If the available space in the bag is smaller than the total going in then divide equally
## print "hitlist = ", hitlist
## print "takestep total_mass = ", total_mass
bag_space = max_load - total_mass
## print "bag_space = ", bag_space
## print "bag_space/float(len(hitlist)) = ", bag_space/float(len(hitlist))
mass_to_add = min(load_add, bag_space/float(len(hitlist)) )
total_mass = total_mass + mass_to_add * float(len(hitlist))
## print "total_mass = ", total_mass, "mass_to_add = ", mass_to_add, "len(hitlist) = ", len(hitlist)
for j in hitlist: #for each outcrop encountered in this step
#print "hit index (j) = ", j
#figure out what type of outcrop it was
for k in hit_ct:
if outcrops[j].id == k.id:
#print "add to type ", k[0]
k.ct += 1
#add to rock bag
k.mass += mass_to_add
def step_dir(ndirs, bdir=vector(0,1,0), biasf=0.0):
"""Determine the direction in which to step
Returns a vector with magnitude of 1.0
bdir = direction of any bias
bias = amount of bias - between 0.0 and 1.0
ndirs = number of possible directions of motion"""
adir = vector(0,1,0)
tmp_dirs = int(ndirs * (1.0 - biasf))
d = float(randint(1-(tmp_dirs)/2,tmp_dirs/2))
# print "bdir a = ", bdir
if bdir <> adir:
#ang =
ang = vector(bdir.x,bdir.y,bdir.z).diff_angle(vector(adir.x, adir.y, adir.z))
# print "ang (degrees) =", ang, ang * 360.0 / 2.0 / pi
newdir = round(ang*ndirs/(2.*pi))
# print "newdir =", ang*ndirs/(2.*pi), "=", newdir
if (adir.x * bdir.y - adir.y * bdir.x) < 0:
bdir = rotate(adir,angle=-(newdir*2.0*pi/float(ndirs)))
else:
bdir = rotate(adir,angle=(newdir*2.0*pi/float(ndirs)))
# print "bdir b = ", bdir
dp = rotate(bdir,angle=d*2.0*pi/float(ndirs))
# if bdir <> adir:
# print "tmp_dirs, d, dp, tmp = ", tmp_dirs, d, dp , ang
return dp
def tri_area(x1, y1, x2, y2, x3, y3):
"""calculate the area of a triangle"""
return x1*y2 - y1*x2 + y1*x3 - x1*y3 + x2*y3 - x3*y2
def stream_crossing(stream_perm,oldpos,newpos,streams,step_size):
"""returns 1 if stream is crossed and move not allowed and 0 otherwise"""
l_cross = 0
## triangle1 = curve(color=(1,0,0))
## triangle1.append(pos=vector(oldpos.x,oldpos.y,0))
## triangle1.append(pos=vector(oldpos.x,oldpos.y,0))
## triangle1.append(pos=vector(oldpos.x,oldpos.y,0))
## triangle1.append(pos=vector(oldpos.x,oldpos.y,0))
## triangle2 = curve(color=(1,0,0))
## triangle2.append(pos=vector(oldpos.x,oldpos.y,0))
## triangle2.append(pos=vector(oldpos.x,oldpos.y,0))
## triangle2.append(pos=vector(oldpos.x,oldpos.y,0))
## triangle2.append(pos=vector(oldpos.x,oldpos.y,0))
if stream_perm < 1.0:
#locate streams close by
#print newpos[0], newpos[1], newpos[2]
for j in streams:
dx = newpos[0] - j.x
dy = newpos[1] - j.y
rad = sqrt(dx * dx + dy * dy) #calculates the distance to each point (outcrop)
##print rad
rad_intercept = less_equal(rad, step_size) #determine if the distance to any point is less than the range of influence (rng)
#rad_list = rad_intercept.tolist()
rad_intercept = rad_intercept * j.crossable #determine if stream is crossable at all
#print rad_intercept
hitlist = sort(nonzero(rad_intercept.flat)).tolist()
for k in hitlist:
x1, y1 = (oldpos.x, oldpos.y)
x2, y2 = (newpos.x, newpos.y)
#look forward along stream
if k < (len(j.x)-1):
x3, y3 = (j.x[k], j.y[k])
x4, y4 = (j.x[k+1], j.y[k+1])
tri1 = tri_area(x1,y1, x2,y2, x3,y3)
tri2 = tri_area(x1,y1, x2,y2, x4,y4)
## triangle1.pos[1] = vector(newpos.x,newpos.y,0)
## triangle1.pos[2] = vector(x3,y3,0)
if tri1/tri2 <= 0.0:
l_cross = 1
break
#look backward along stream
if k > 0:
x3, y3 = (j.x[k-1], j.y[k-1])
x4, y4 = (j.x[k], j.y[k])
tri1 = tri_area(x1,y1, x2,y2, x3,y3)
tri2 = tri_area(x1,y1, x2,y2, x4,y4)
if tri1/tri2 <= 0.0:
l_cross = 1
break
if l_cross == 1:
break
if l_cross == 1: #then check to see if you can cross stream
if stream_perm < random():
return 1
else:
return 0
else:
return 0
#end streams
def topo_bias(nslope,sslope,eslope,wslope, xmin, ymax, dx, pos, dp, max_slope, l_slope = 0.0):
if l_slope <> 0.0:
#determine cell location
r = int(((ymax - pos.y) / dx) - (dx/2.0))
c = int(((pos.x - xmin) / dx) - (dx/2.0))
#determine applicable slope
dir = norm(dp)
if (-0.5 < dir.x <= 0.5): #then move is north or south
if dir.y >= 0: #move is north
s = nslope.data[r,c]
else: #move is south
s = sslope.data[r,c]
else: #move is east or west
if dir.x > 0: #move is east
s = eslope.data[r,c]
else: #move is west
s = wslope.data[r,c]
#evaluate probability of not making the move
p_not = (abs(s) / max_slope) * l_slope # should give a value less than 1; greater slope = greater reluctance
if random() < p_not:
#print "don't do it"
return 1
else:
return 0
else:
return 0
def attract_grad(pos, oray):
#oray is 2d array with id, x, y, Azero, halflife: for locations of outcrops
#r, c = arr.shape
dx = pos[0] - oray[:,1]
dy = pos[1] - oray[:,2]
rad = sqrt(dx * dx + dy * dy) #calculates the distance to each point (outcrop)
gexp = -0.693 * rad / oray[:,5]
gray = oray[:,6] * ((-0.693 * oray[:,4] / oray[:,5]) * exp(gexp))
norm_x = dx / (dx + dy)
norm_y = dy / (dx + dy)
grayx = sum(gray * norm_x)
grayy = sum(gray * norm_y)
#print "attractor gradient =", vector(grayx, grayy, 0.0)
return vector(grayx, grayy, 0.0)
##def write_walk(hit_ct, outfile):
## outln = ""
## for i in hit_ct:
## outln = outln + str(i.ct) + " "
## outln = outln + str(i.mass) + " "
## outln = outln + "\n"
## outfile.write(outln)
#attract = []
cherttypes = []
#####################################################
# INPUT PARAMETERS
#####################################################
#RANGE (rng) IS THE RANGE OF INFLUENCE OF EACH LOCATION
rng = 250
#NUMBER OF STEPS IN RANDOM WALK (nsteps)
nsteps = 100
step_size = 1000
ndirs = 360 #the number of directions you can go off in
#NUMBER OF WALKS (nwalks)
nwalks = 20
#RANDOM NUMBER GENERATION SPECIFICATIONS
seed = 1
#MAXIMUM LOAD CARRIED
max_load = 20.
min_load = 0.0
#mass picked up from each outcrop encountered
load_add = 1.0
#load attrition rate (as a fraction) per rocktype per step
load_atr = 0.01
#test for convergence using an overnight run set l_converg(start,stop)
# the simulation will repeat l_converg times
# the maximum number of walks will be nwalks * 2^l_converg
# To have only one set of walks run set
# l_converg = range(0,1)
# l_scale is the base to which the number of simulation is raised should be an integer
# eg will give nwalks*l_scale**l_converg[i]
l_converg = range(0,1)
l_scale = 1
#draw 3d output if l_draw = 1
l_draw = 1
# draw histogram if l_draw is also 1
l_draw_graph = 1
# draw dem
l_draw_topo = 1
#scale elevation data
z_scale = 10.0 # to scale the elevation data (0.0 means a 2d map)
#setup scene window
if l_draw == 1:
sw = 500
sh = 800
sx = 0
sy = 0
scene.width = sw
scene.height = sh
scene.x = sx
scene.y = sy
#site locations
sites = []
col = vector(1,1,0)
sites.append(outcrop("la-cote", 639495.5157, 3311775.908, l_draw, 2.*rng, col))
sites.append(outcrop("Le Moustier", 682537.139, 3298763.248, l_draw, 2.*rng, col))
sites.append(outcrop("Le Flag", 683737.1246, 3282782.176, l_draw, 2.*rng, col))
sites.append(outcrop("Caminade", 697025.6269, 3283652.463, l_draw, 2.*rng, col))
#site where random walks originate - must be exactly the same as one of the location entered above (in the list sites[])
start_site = "Le Moustier"
# the home_sickness value restricts motion opposite the direction of home
# 0.0 <= home_sickness <= 1.0
# if home_sickness = 0.5 then the chert-man can only move in a direction that is within 90 degrees of the home direction.
home_sickness = 0.5 #bias toward heading home after the walk (must be between 0 and 1.0)
# stream permeability -
# 0.0 <= stream_perm <= 1.0
# stream_perm is the probability of being able to cross a stream
stream_perm = 1.0
# to consider slope in the program set l_slope to a value greater than 0.0
# if l_slope > 1.0 then that sets an upper limit as to the slope
# that chertman can scale.
l_slope = 0.0
#l_stuck is the number of times for chertman to try and fail to escape somewhere.
l_stuck = 100
# set up attractors for the different directions of motion
# l_attract = 0 or 1
l_attract = 0
# you need to know ahead of time the ID's of the data
# to enter attractor('id', 'Azero', 'halflife', 'attractive', 'atr_rate', 'prob_pickup')
# prob_pickup is a value between 0.0 and 1.0: Default = 1.0
# it tells the probability that a type of chert will be picked up
# if encountered
cherttypes.append(cherttype('1', Azero = 100.0, halflife = 1000, attractive = 1.0,
atr_rate = load_atr, prob_pickup = 1.0))
cherttypes.append(cherttype('2', atr_rate = load_atr, prob_pickup = 1.0))
cherttypes.append(cherttype('3', atr_rate = load_atr, prob_pickup = 1.0))
cherttypes.append(cherttype('4', atr_rate = load_atr, prob_pickup = 1.0))
cherttypes.append(cherttype('5', atr_rate = load_atr, prob_pickup = 1.0))
cherttypes.append(cherttype('7', atr_rate = load_atr, prob_pickup = 1.0))
cherttypes.append(cherttype('8', atr_rate = load_atr, prob_pickup = 1.0))
cherttypes.append(cherttype('25', atr_rate = load_atr, prob_pickup = 1.0))
norm_abias = 1.0
# INPUT FILES
# Outcrop data as ordered list of data points with header line and format
# X Y ID
f_outcrop = "rechecked_xy_transf.txt"
# Stream data as an ordered list of data points with headerline and format
# X Y StreamName
f_stream = "streams_transformed2.txt"
# topography
f_topo = "rastert_asciito1.txt"
#####################################################
# END OF INPUT PARAMETERS
#####################################################
print "File Name = ", argv[0]
wfout_file = basename(argv[0]) + ".dat"
#output filenames
#full data
#wfout_file =
for i in sites:
if i.id == start_site:
start_pos = vector(i.x, i.y, 0.)
print i.id, i.x, i.y
#set scene center
if l_draw == 1:
scene.center = start_pos
scene.background = (1,1,1)
scene.autoscale = 0
scene.range = 50000
if l_draw_graph == 1:
gdisp = gdisplay(x=sw, y=0, width=sw, height=sh/2)
## ytitle="Average Mass (every 10 steps)", xtitle="walk",
## ymin=min_load, ymax=max_load)
masshist = ghistogram(bins=arange(min_load, max_load, load_add), color=color.red)
# accumulate=1, average=1)
## Tcurve = gcurve(color=color.cyan) # a connected curve object
## Tdots = gdots(color=(1,1,0))
''' READ IN OUTCROP DATA '''
outcrops = [] #list with locations of outcrops and their type
oray = [] #array for outcrops giving columns ['id','x','y','z','Azero','halflife','attractive'
# First open the file to read(r)
infile = open(f_outcrop,"r")
# read the file into a list then print each item
ncrops = 0
crop_ct = 0
hit_ct = []
attract_array = [] # the first column will be the attractor depth and the second the half life
old_crop = ""
filelines = infile.readlines()
for i in range(1, len(filelines)): #skip first header line
ncrops = ncrops+1
l = split(strip(filelines[i]))
dat_id = l[2]
dat_x = l[0]
dat_y = l[1]
#print "dat_id, old_crop", dat_id, old_crop
#print "filelines = ", filelines[i]
same_as = 0
for j in outcrops:
#print "i, j, j[0], old_crop = ", i, j, j[0], old_crop
if j.id == dat_id: #j[0] is the outcrop id
same_as = 1
#print same_as
if same_as == 0:
crop_ct += 1
old_crop = dat_id
#print "old_crop, line =", old_crop
hit_ct.append(hit_counter(dat_id))
outcrops.append(outcrop(dat_id, dat_x, dat_y, l_draw, rng))
#oray.append((float(dat_id), float(dat_x), float(dat_y), 0.0, 0.0))
if l_attract <> 0:
for j in cherttypes:
if dat_id == j.id:
#print "oray[-1] = ", oray[-1][3]
#print "j.Azero = ", j.Azero
#print "j.halflife = ", j.halflife
oray.append((float(dat_id), float(dat_x), float(dat_y), 0.0, j.Azero, j.halflife, j.attractive))
#oray[-1][3] = 10.0
#oray[-1][4] = j.halflife
else:
oray.append((float(dat_id), float(dat_x), float(dat_y)))
## if l_attract == 1:
## if len(l) >= 6:
## outcrops[-1].Azero = float(l[5])
## attract_array.append(outcrops[-1].Azero)
## if len(l) >= 7:
## outcrops[-1].halflife = float(l[6])
## attract_array.append(outcrops[-1].halflife)
for j in cherttypes:
if j.id == dat_id:
j.x.append(float(dat_x))
j.y.append(float(dat_y))
# Now close it again
infile.close()
oray = array(oray)
for i in hit_ct:
for j in cherttypes:
if i.id == j.id:
i.atr_rate = j.atr_rate
i.prob_pickup = j.prob_pickup
#set outcrop color
for i in outcrops:
if i.id == "1":
i.ball.color = vector(0.0,0.0,0.0)
if i.id == "2":
i.ball.color = vector(.83,.6,.4)
if i.id == "3":
i.ball.color = vector(.7,.0,.4)
if i.id == "4":
i.ball.color = vector(0.0,0.8,0.1)
if i.id == "5":
i.ball.color = vector(0.9,0.9,.15)
if i.id == "25":
i.ball.color = vector(0.99,0.6,.0)
if i.id == "8":
i.ball.color = vector(.7,.52,.4)
print "Number of outcrops = ", ncrops
#print "oray "
#print oray
#print "oray first row"
#print oray[0,:]
##print "oc"
##print oc
print "number of different types of outcrops = ", crop_ct
#print hit_ct
'''STREAMS'''
streams = [] #list of streams
oc = []
# First open the file to read(r)
infile = open(f_stream,"r")
# read the file into a list then print each item
ct = 0
str = -1
str_old = ""
flines = infile.readlines()
for i in range(1, len(flines)):
l = split(strip(flines[i]))
dat_id = l[2]
dat_x = l[0]
dat_y = l[1]
dat_size = l[3]
if str_old <> dat_id:
str_old =dat_id
streams.append(stream(dat_id))
str = str+1
#print streams[str].name
streams[str].x.append(float(dat_x))
streams[str].y.append(float(dat_y))
streams[str].z.append(0.0)
streams[str].size.append(int(dat_size))
# Now close file
infile.close()
for i in streams:
i.x = array(i.x)
i.y = array(i.y)
i.z = array(i.z)
i.size = array(i.size)
i.cross()
#rotate streams, outcrops, sites
origin = vector(680097, 3294948, 0)
rangle = 2.0
for i in streams:
for j in range(len(i.x)):
#print i.name, i.x[j], i.y[j]
posvec = vector(i.x[j], i.y[j], 0.0)
posvec = posvec - origin #translate
posvec = rotate(posvec, angle=rangle*pi/180, axis=(0,0,1)) #rotate
posvec = posvec + origin #translate
i.x[j] = posvec.x
i.y[j] = posvec.y
#print " ", i.x[j], i.y[j]
for i in outcrops:
posvec = vector(float(i.x), float(i.y), 0.0)
posvec = posvec - origin #translate
posvec = rotate(posvec, angle=rangle*pi/180, axis=(0,0,1)) #rotate
posvec = posvec + origin #translate
i.x = posvec.x
i.y = posvec.y
for i in sites:
posvec = vector(float(i.x), float(i.y), 0.0)
posvec = posvec - origin #translate
posvec = rotate(posvec, angle=rangle*pi/180, axis=(0,0,1)) #rotate
posvec = posvec + origin #translate
i.x = posvec.x
i.y = posvec.y
'''TOPOGRAPHY'''
#READ IN TOPOGRAPHY FROM arcinfo GIS OUTPUT OF RASTER FILE (ASCII)
### reads in the header information also
topo = raster_import(f_topo)
print "TOPOGRAPHY"
print "nrows, ncols, dx, xmin, ymax"
print topo.nrows, topo.ncols, topo.dx, topo.xmin, topo.ymax
print topo.data[0:5,0:20]
#view(topo)
infile.close()
#Elevate outcrop data
for i in cherttypes:
## print "cherttype = ", i.id
## print "x in"
## print i.x
for j, k in zip(i.x, i.y):
elev = elevate(vector(j,k,0.0),topo)
i.z.append(elev.z)
i.makearrays()
## print "x out"
## print i.x
## print "y"
## print i.y
## print "z"
## print i.z
#calculate slopes in 4 directions
if l_slope > 0.0:
tmp = ( topo.data[0:nrows-1,:] - topo.data[1:nrows,:] ) / dx
nslope = zeros((nrows*ncols,), Float)
nslope = resize(nslope, (nrows, ncols))
nslope[1:nrows,:] = tmp
## for i in range(1,nrows):
## for j in range(ncols):
## nslope[i,j] = tmp[i-1,j]
tmp = ( topo.data[1:nrows,:] - topo.data[0:nrows-1,:] ) / dx
sslope = zeros((nrows*ncols,), Float)
sslope = resize(sslope, (nrows, ncols))
sslope[0:nrows-1,:] = tmp
## for i in range(0,nrows-1):
## for j in range(ncols):
## sslope[i,j] = tmp[i,j]
tmp = ( topo.data[:,1:ncols] - topo.data[:,0:ncols-1] ) / dx
eslope = zeros((nrows*ncols,), Float)
eslope = resize(eslope, (nrows, ncols))
eslope[:,0:ncols-1] = tmp
## for i in range(nrows):
## for j in range(0,ncols-1):
## eslope[i,j] = tmp[i,j]
tmp = ( topo.data[:,0:ncols-1] - topo.data[:,1:ncols] ) / dx
wslope = zeros((nrows*ncols,), Float)
wslope = resize(wslope, (nrows, ncols))
wslope[:,1:ncols] = tmp
## for i in range(nrows):
## for j in range(0,ncols-1):
## wslope[i,j] = tmp[i,j-1]
max_slope = max(max(nslope))
nslope = raster_map(nslope, xmin, ymax, dx, nrows, ncols)
sslope = raster_map(sslope, xmin, ymax, dx, nrows, ncols)
eslope = raster_map(eslope, xmin, ymax, dx, nrows, ncols)
wslope = raster_map(wslope, xmin, ymax, dx, nrows, ncols)
#nslope.draw_map( 120, 200, 240, 340)
print "MAXIMUM SLOPE IN DATASET (max_slope) = ", max_slope
print "Maximum slope chert-man can scale = ", max_slope / l_slope
##print "wSLOPE"
##print wslope[0:5,0:20]
''' CALCULATE ATTRACTOR SLOPE MORPHOLOGY '''
# attract_array has columns [Azero, halflife, active]
# construct raster array coincident with the topography
if l_attract <> 0:
print 'calculating attractor slope morphology'
attx = resize(zeros((nrows*ncols,), Float), (nrows,ncols))
atty = resize(zeros((nrows*ncols,), Float), (nrows,ncols))
zp = 0.0
for i in range(nrows):
print "row = %i of %i" % (i, nrows)
yp = ymax - dx*i
for j in range(ncols):
xp = xmin + dx*j
#print "xp, yp, zp = ", xp, yp, zp
attx[i,j] = attract_grad(vector(xp,yp,0.0),oray).x
atty[i,j] = attract_grad(vector(xp,yp,0.0),oray).y
## print "attract_array"
## print attract_array
# create raster map for attractors
attx_map = raster_map(attx, xmin, ymax, dx, nrows, ncols)
atty_map = raster_map(atty, xmin, ymax, dx, nrows, ncols)
#PLOT TO VPYTHON DISPLAY
if l_draw == 1:
print "Drawing ..."
## #plot sites
## for i in sites:
## print i.x, i.y
## sphere(pos=(i.x, i.y), radius=rng*2.0, color=(1,1,0))
##
## #plot outcrops
## col = vector(1,0,1)
## for i in outcrops:
## sphere(pos=(i.x, i.y), radius=rng, color=col)
#start_pos = elevate(start_pos, topo, z_scale)
#plot outcrop elevations
for i in sites:
i.elevate(topo, z_scale)
for i in outcrops:
i.elevate(topo, z_scale)
#plot streams
print "STREAM NAME Number of points in stream"
for i in streams:
i.elevate(topo, z_scale)
print " ", i.name, len(i.x)
river = curve(color = (0,0,1))
for j in range(len(i.x)):
river.append(pos = (i.x[j], i.y[j], i.z[j]))
#tx = str(j)
#label each point on line
## label(pos=(float(i.x[j]), float(i.y[j])), text=repr(j),
## box=0, opacity = 0, xoffset=-5, yoffset=5, height=10, color=color.blue)
#print float(i.x[j]), float(i.y[j])
#rate(1)
if l_draw_topo == 1:
#plot topography
# resize topo array to area of interest
#topo.draw_map( 60, 290, 200, 380)
#topo.draw_map( 100, 200, 250, 350, scale = z_scale)
#topo.draw_map(scale = z_scale)
topo.line_3d(100, 200, 250, 350,scale = z_scale)
if l_attract <> 0:
attx_map.line_3d( 60, 290, 200, 380)
'''CALCULATE DISTRIBUTIONS BASED ON DISTANCE OF SITE FROM OUTCROPS'''
# Distance weighted
dist_fout = open("radius_based_distr.txt","w")
dist_fout.write("Calculating distance weighted distributions \n")
for i in sites:
outline = "Site =,%s" % i.id + "\n"
dist_fout.write(outline)
type_dist = []
inv_rad = []
inv_rad2 = []
for j in cherttypes:
dx = j.x - i.x
dy = j.y - i.y
rad = sqrt(dx * dx + dy * dy)
wt = add.reduce(rad)
type_dist.append(wt)
irad = add.reduce(1.0/rad)
inv_rad.append(irad)
irad2 = add.reduce(1.0/(rad*rad))
inv_rad2.append(irad2)
tot_dist = reduce(add, type_dist)
tot_inv = reduce(add, inv_rad)
tot_inv2 = reduce(add, inv_rad2)
outline = "Chert type, Total distance to type, Distance weight, Inverse distance Weight, Inverse Distance Squared \n"
dist_fout.write(outline)
type_frac = []
inv_frac = []
inv2_frac = []
print "id, total radius, wt, total inverse distance, wt"
for j, k in enumerate(cherttypes):
type_frac.append(type_dist[j] / tot_dist)
inv_frac.append(inv_rad[j] / tot_inv)
inv2_frac.append(inv_rad2[j] / tot_inv2)
outline = "%s, %10.20f, %10.20f, %10.20f, %10.20f" % (k.id, type_dist[j], type_frac[-1], inv_frac[-1], inv2_frac[-1]) + "\n"
dist_fout.write(outline)
print k.id, type_dist[j], type_frac[-1], inv_rad[j], inv_frac[-1], inv2_frac[-1]
tot_frac = reduce(add, type_frac)
tot_invfrac = reduce(add, inv_frac)
tot_invfrac2 = reduce(add, inv2_frac)
outline = "Total, %10.20f, %1.5f, %1.5f, %1.5f" % (tot_dist, tot_frac, tot_invfrac, tot_invfrac2) + "\n"
dist_fout.write(outline)
dist_fout.close()
##for i in sites:
## dx = oray[:,1] - i.x
## dy = oray[:,2] - i.y
## rad = sqrt(dx * dx + dy * dy)
## tot_dist = add.reduce(rad)
## print "tot_dist = ", tot_dist
## for j in cherttypes:
## type_rad = 0.0
## for k in outcrops:
## if j.id == k.id:
## dx = i.x - k.x
## dy = i.y - k.y
## type_rad = type_rad + (dx * dx + dy * dy) ** 0.5
##
##
#output data for every walk
wfout = open(wfout_file,"a")
outline = "RANDOM WALK SIMULATION - " + datetime.today().strftime("%A %B %d, %Y. - %I:%M:%S") + "\n"
wfout.write(outline)
outline = "File Name = %s" % argv[0] + "\n"
wfout.write(outline)
outline = "Starting Site (start_site) = %s" % start_site + "\n"
wfout.write(outline)
outline = "RANGE OF INFLUENCE OF EACH OUTCROP (rng) = %1.1f" % rng + "\n"
wfout.write(outline)
outline = "NUMBER OF STEPS IN RANDOM WALK (nsteps) = %i" % nsteps + "\n"
wfout.write(outline)
outline = "SIZE OF EACH RANDOM STEP (step_size) = %1.1f" % step_size + "\n"
wfout.write(outline)
outline = "NUMBER OF MOVEMENT DIRECTIONS (ndirs) = %i" % ndirs + "\n"
wfout.write(outline)
outline = "NUMBER OF WALKS IN SIMULATION (nwalks) = %i" % nwalks + "\n"
wfout.write(outline)
outline = "MAXIMUM LOAD (max_load) = %1.1f kg" % max_load + "\n"
wfout.write(outline)
outline = "MINIMUM LOAD (min_load) = %1.1f kg" % min_load + "\n"
wfout.write(outline)
outline = "MASS PICKED UP per step (load_add) = %1.1f kg" % load_add + "\n"
wfout.write(outline)
outline = "ATTRITION RATE per step (load_atr) = \n"
wfout.write(outline)
outline = ''
for i in hit_ct:
outline = outline + "ID=%s atr_rate=%1.3f ;" % (i.id, i.atr_rate)
outline = outline + "\n"
wfout.write(outline)
outline = "Probability of crossing stream (stream_perm)= %1.3f " % stream_perm + "\n"
wfout.write(outline)
outline = "Slope effects (l_slope) = %1.3f " % l_slope + "\n"
wfout.write(outline)
outline = "Maximum escape attempts if stuck (l_stuck) = %i " % l_stuck + "\n"
wfout.write(outline)
outline = ""
for i in hit_ct:
tmp1 = "Nhits-ID=%s " % i.id #the type of outcrops encountered
outline = outline + tmp1 #+ tmp2
for i in hit_ct:
tmp1 = "Mass-ID=%s " % i.id #the type of outcrops encountered
outline = outline + tmp1 #+ tmp2
outline = outline + "walk_steps "
outline = outline + "Total\n"
wfout.write(outline)
chertman = sphere(radius = rng, color=color.red, pos=start_pos)
chertman.pos = elevate(chertman.pos, topo, z_scale)
#pathway = curve(color = (0,1,0))
#random walk
tot_mass = [] #list giving the total mass at the end of each walk
newpos = vector(0,0,0)
tot_walks = 0
adir = vector(0,1,0)
abias = 0.0
nwalks_init = nwalks
for lc in l_converg: # for convergence testing
nwalks = nwalks_init * (l_scale**lc)
print "nwalks = ", nwalks
nhits = 0
for w in range(1,nwalks+1):
tot_walks += 1
# reset bag contents to zero for each walk
for i in hit_ct:
i.mass = 0.0
total_mass = 0.0
if l_draw == 1:
#pathway.visible = 0
pathway = curve(color = (0,1,0), radius=100)
pathway.append(pos=elevate(start_pos, topo, z_scale))
route = []
route.append(start_pos)
## print "start"
## rate(.1)
## print "go"
'''HEAD OUT'''
walk_steps = 0
for i in range(1,nsteps+1):
#attrition of material over step
total_mass = 0.0
for i in hit_ct:
#print "item mass = ", i.id, i.mass
i.mass = i.mass - (i.mass * i.atr_rate) #mass attrition function
i.mass = max(i.mass, 0.0) #can't have negative mass
total_mass = total_mass + i.mass
l_stuck_test = 0
l_again = 1
while l_again == 1:
l_again =0
#determine new position
# adjust for attractors by determining bias direction
if l_attract <> 0:
agrad = attract_grad(route[-1], oray)
if agrad.x <> 0 or agrad.y <> 0:
adir = norm(agrad)
abias = min(mag(agrad)/norm_abias, 1.0)
else:
abias = 0.0
print "adir, abias = ", adir, abias
#take step
dp = step_size * step_dir(ndirs, adir, abias)
newpos = route[-1]+dp
chertman.pos = elevate(newpos, topo, z_scale)
#adjust for streams
l_again = stream_crossing(stream_perm,route[-1],newpos,streams,step_size)
#adjust for topography (slope)
if l_again == 0 and l_slope > 0.0:
l_again = topo_bias(nslope,sslope,eslope,wslope,
xmin, ymax, dx, route[-1], dp, max_slope, l_slope)
#see that chertman is not stuck in one place forever
if l_again ==1:
l_stuck_test += 1
if l_stuck_test >= l_stuck:
l_again = 0
walk_steps += 1
if l_draw == 1:
pathway.append(pos=elevate(newpos, topo, z_scale))
## rate(.2)
route.append(newpos)
# Take step and add appropriate mass
#print "o total_mass = ", total_mass
take_step(nhits,ndirs,step_size,oray,hit_ct,outcrops,load_add, max_load, total_mass)
'''HEAD HOME'''
if l_draw == 1:
#pathway.append(color = (1,0,0))
pathway.color=(1,0,0)
ct = nsteps
while mag(newpos - start_pos) > step_size:
#attrition of material over step
total_mass = 0.0
for i in hit_ct:
#print "item mass = ", i.id, i.mass
i.mass = i.mass - (i.mass * i.atr_rate) #mass attrition function
i.mass = max(i.mass, 0.0) #can't have negative mass
total_mass = total_mass + i.mass
l_stuck_test = 0
l_again = 1
while l_again == 1:
l_again = 0
#determine new position
home_dir = start_pos - newpos
home_dir = norm(vector(home_dir.x, home_dir.y, 0.0))
## print "start_pos = ", start_pos
## print "newpos = ", newpos
## print " ct = ", ct
## print "home_dir = ", home_dir
dp = step_size * step_dir(ndirs, home_dir, home_sickness) #move biased by 0.5 in the direction of home
newpos = route[-1]+dp
chertman.pos = elevate(newpos, topo, z_scale)
#adjust for streams
l_again = stream_crossing(stream_perm,route[-1],newpos,streams,step_size)
#adjust for topography (slope)
if l_again == 0 and l_slope > 0.0:
l_again = topo_bias(nslope,sslope,eslope,wslope,
xmin, ymax, dx, route[-1], dp, max_slope, l_slope)
#see that chertman is not stuck in one place forever
if l_again ==1:
l_stuck_test += 1
if l_stuck_test >= l_stuck:
l_again = 0
walk_steps += 1
if l_draw == 1:
pathway.append(pos=elevate(newpos, topo, z_scale))
route.append(newpos)
#print "i total_mass = ", total_mass
take_step(nhits,ndirs,step_size,oray,hit_ct,outcrops,load_add, max_load, total_mass)
ct += 1
if l_draw == 1:
pathway.append(pos=elevate(start_pos, topo, z_scale))
# print "pause 2"
#record mass at end of walk
total_mass = 0.0
for i in hit_ct:
tmp_mass = i.avmass * (tot_walks - 1) + i.mass
i.avmass = tmp_mass / float(tot_walks)
i.stmass.append(i.mass)
total_mass = total_mass + i.mass
#print " id = ", i.id, ": mass = ", i.mass
tot_mass.append(total_mass)
#write out at end of each walk
outline = ""
for i in hit_ct:
outline = outline + "%i " % i.ct
for i in hit_ct:
outline = outline + "%1.5f " % i.mass
outline = outline + "%i " % walk_steps
outline += "\n"
wfout.write(outline)
if (w % 10) == 0:
print "Walk ", w, " / ", nwalks, ": tot_mass[-10:-1] = ", tot_mass[-10:-1]
#print "Walk =", tot_walks, " total_mass = ", tot_mass[-1]
if l_draw_graph == 1:
masshist.plot(data=tot_mass) # plot the Total Mass distribution
## Tcurve.plot(pos=(tot_walks, mean(tot_mass)))
## Tdots.plot(pos=(tot_walks, mean(tot_mass[-10:-1])))
print "number of steps = ", nsteps
print "number of hits = ", nhits
print "types of hits"
wfout.close()
for i in hit_ct:
#calculate statistics
if len(i.stmass) > 1:
print "i.stmass = ", i.stmass
i.stdev = standardDeviation(i.stmass)
else:
i.stdev = 0.0
i.min = min(i.stmass)
i.max = max(i.stmass)
i.median = median(i.stmass)
print "id = ", i.id, ": nhits = ", i.ct, ": Average mass = ", i.avmass
#i.skew = skewness(i.stmass)
print "Average total mass at end of run = ", mean(tot_mass)
#Write output file
fout = open("random_walk_out.txt","a")
fout.write("\n")
outline = "RANDOM WALK SIMULATION COMPLETED - " + datetime.today().strftime("%A %B %d, %Y. - %I:%M:%S") + "\n"
fout.write(outline)
outline = "File Name = %s" % argv[0] + "\n"
fout.write(outline)
outline = "Starting Site (start_site) = %s" % start_site + "\n"
fout.write(outline)
outline = "RANGE OF INFLUENCE OF EACH OUTCROP (rng) = %1.1f" % rng + "\n"
fout.write(outline)
outline = "NUMBER OF STEPS IN RANDOM WALK (nsteps) = %i" % nsteps + "\n"
fout.write(outline)
outline = "SIZE OF EACH RANDOM STEP (step_size) = %1.1f" % step_size + "\n"
fout.write(outline)
outline = "NUMBER OF MOVEMENT DIRECTIONS (ndirs) = %i" % ndirs + "\n"
fout.write(outline)
outline = "NUMBER OF WALKS IN SIMULATION (nwalks) = %i" % nwalks + "\n"
fout.write(outline)
outline = "MAXIMUM LOAD (max_load) = %1.1f kg" % max_load + "\n"
fout.write(outline)
outline = "MINIMUM LOAD (min_load) = %1.1f kg" % min_load + "\n"
fout.write(outline)
outline = "MASS PICKED UP per step (load_add) = %1.1f kg" % load_add + "\n"
fout.write(outline)
outline = "ATTRITION RATE per step (load_atr) = %1.3f " % load_atr + "\n"
fout.write(outline)
outline = ''
for i in hit_ct:
outline = outline + "ID=%s atr_rate=%1.3f ;" % (i.id, i.atr_rate)
outline = outline + "\n"
fout.write(outline)
outline = "Probability of crossing stream (stream_perm)= %1.3f " % stream_perm + "\n"
fout.write(outline)
outline = "Slope effects (l_slope) = %1.3f " % l_slope + "\n"
fout.write(outline)
outline = "Maximum escape attempts if stuck (l_stuck) = %i " % l_stuck + "\n"
fout.write(outline)
outline = "Average total mass (mean(tot_mass)) = %1.3f kg" % mean(tot_mass) + "\n"
fout.write(outline)
fout.write("\n")
outline = " ID# #HITS Mass Average_mass Standard_Deviation Minimum Maximum Median \n"
fout.write(outline)
for i in hit_ct:
outline = " %5s %5i %5.5f %5.5f %5.5f %5.5f %5.5f %5.5f" % (i.id, i.ct, i.mass, i.avmass, i.stdev, i.min, i.max, i.median) + "\n"
fout.write(outline)
if len(tot_mass) > 1:
totstdev = standardDeviation(tot_mass)
else:
totstdev = 0.0
outline = " %5s %5i %5.5f %5.5f %5.5f %5.5f %5.5f %5.5f" % ("Total", 0, 0, mean(tot_mass),totstdev, min(tot_mass), max(tot_mass), median(tot_mass)) + "\n"
fout.write(outline)
fout.close()
## print rad_list
if (len(argv) > 1):
os._exit(0)
