Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
mg = fr.route_flow()
# mg.calculate_gradient_across_cell_faces(mg.at_node['topographic__elevation'])
#neighbor_slopes = mg.calculate_gradient_along_node_links(mg.at_node['topographic__elevation'])
#mean_slope = np.mean(np.fabs(neighbor_slopes),axis=1)
#max_slope = np.max(np.fabs(neighbor_slopes),axis=1)
#mg,_,capacity_out = tl.erode(mg,dt,slopes_at_nodes='topographic__steepest_slope')
#mg,_,capacity_out = tl.erode(mg,dt,slopes_at_nodes=max_slope)
mg_copy = deepcopy(mg)
mg, _ = tl.erode(mg, dt, stability_condition='loose')
if i % 20 == 0:
print('loop ', i)
print('subdivisions of dt used: ', tl.iterations_in_dt)
print('max_slope', np.amax(
mg.at_node['topographic__steepest_slope'][mg.core_nodes]))
pylab.figure("long_profiles")
profile_IDs = prf.channel_nodes(mg, mg.at_node['topographic__steepest_slope'],
mg.at_node['drainage_area'], mg.at_node['flow_receiver'])
dists_upstr = prf.get_distances_upstream(mg, len(mg.at_node['topographic__steepest_slope']),
profile_IDs, mg.at_node['links_to_flow_receiver'])
prf.plot_profiles(dists_upstr, profile_IDs, mg.at_node[
'topographic__elevation'])
# mg.update_boundary_nodes()
#vid.add_frame(mg, 'topographic__elevation')
print('Completed the simulation. Plotting...')
time_off = time()
#Finalize and plot
elev = mg['node']['topographic__elevation']
#neighbor_slopes = mg.calculate_gradient_along_node_links(mg.at_node['topographic__elevation'])
#mean_slope = np.mean(np.fabs(neighbor_slopes),axis=1)
#max_slope = np.max(np.fabs(neighbor_slopes),axis=1)
#mg,_,capacity_out = tl.erode(mg,dt,slopes_at_nodes='topographic__steepest_slope')
#mg,_,capacity_out = tl.erode(mg,dt,slopes_at_nodes=max_slope)
mg_copy = deepcopy(mg)
mg, _ = tl.erode(mg, dt, stability_condition='loose')
if i % 20 == 0:
print('loop ', i)
print('subdivisions of dt used: ', tl.iterations_in_dt)
print('max_slope', np.amax(
mg.at_node['topographic__steepest_slope'][mg.core_nodes]))
pylab.figure("long_profiles")
profile_IDs = prf.channel_nodes(mg, mg.at_node['topographic__steepest_slope'],
mg.at_node['drainage_area'], mg.at_node['flow_receiver'])
dists_upstr = prf.get_distances_upstream(mg, len(mg.at_node['topographic__steepest_slope']),
profile_IDs, mg.at_node['links_to_flow_receiver'])
prf.plot_profiles(dists_upstr, profile_IDs, mg.at_node[
'topographic__elevation'])
# mg.update_boundary_nodes()
#vid.add_frame(mg, 'topographic__elevation')
print('Completed the simulation. Plotting...')
time_off = time()
#Finalize and plot
elev = mg['node']['topographic__elevation']
#imshow.imshow_node_grid(mg, elev)
#mg,_,capacity_out = tl.erode(mg,dt,slopes_at_nodes='topographic__steepest_slope')
#mg,_,capacity_out = tl.erode(mg,dt,slopes_at_nodes=max_slope)
mg_copy = deepcopy(mg)
mg, _ = sde.erode(mg, dt)
# print sde.iterations_in_dt
# print 'capacity ', np.amax(capacity_out[mg.core_nodes])
# print 'rel sed ',
# np.nanmax(sed_in[mg.core_nodes]/capacity_out[mg.core_nodes])
if i % 100 == 0:
print('loop ', i)
print('max_slope', np.amax(
mg.at_node['topographic__steepest_slope'][mg.core_nodes]))
pylab.figure("long_profiles")
profile_IDs = prf.channel_nodes(mg, mg.at_node['topographic__steepest_slope'],
mg.at_node['drainage_area'], mg.at_node['flow__receiver_node'])
dists_upstr = prf.get_distances_upstream(mg, len(mg.at_node['topographic__steepest_slope']),
profile_IDs, mg.at_node['flow__link_to_receiver_node'])
prf.plot_profiles(dists_upstr, profile_IDs, mg.at_node[
'topographic__elevation'])
#vid.add_frame(mg, 'topographic__elevation')
print('Completed the simulation. Plotting...')
time_off = time()
#Finalize and plot
elev = mg['node']['topographic__elevation']
#imshow.imshow_node_grid(mg, elev)
print('Done.')
print("Short step!")
dt = time_to_run - elapsed_time
mg = fr.run_one_step()
#print 'Area: ', numpy.max(mg.at_node['drainage_area'])
#mg = fsp.erode(mg)
mg,_,_ = sp.erode(mg, dt, K_if_used='K_values')
#plot long profiles along channels
if numpy.allclose(elapsed_time%out_tstep,0.) or numpy.allclose(elapsed_time%out_tstep,1.):
pylab.figure("long_profiles")
profile_IDs = prf.channel_nodes(mg, mg.at_node['topographic__steepest_slope'],
mg.at_node['drainage_area'], mg.at_node['flow__receiver_node'])
dists_upstr = prf.get_distances_upstream(mg, len(mg.at_node['topographic__steepest_slope']),
profile_IDs, mg.at_node['flow__link_to_receiver_node'])
prf.plot_profiles(dists_upstr, profile_IDs, mg.at_node['topographic_elevation'])
x_profiles.append(dists_upstr[:])
z_profiles.append(mg.at_node['topographic_elevation'][profile_IDs])
S_profiles.append(mg.at_node['steepest_slope'][profile_IDs])
A_profiles.append(mg.at_node['drainage_area'][profile_IDs])
prf.plot_profiles(dists_upstr, profile_IDs, mg.at_node['topographic__elevation'])
#add uplift
mg.at_node['topographic__elevation'][mg.core_nodes] += 5.*uplift*dt
elapsed_time += dt
#Finalize and plot
elev = mg['node']['topographic__elevation']
elev_r = mg.node_vector_to_raster(elev)
dt=10000.
out_tstep = 50000.
elapsed_time = 0. #total time in simulation
while elapsed_time < time_to_run:
if elapsed_time+dt>time_to_run:
print("Short step!")
dt = time_to_run - elapsed_time
mg = fr.run_one_step()
#print 'Area: ', numpy.max(mg.at_node['drainage_area'])
#mg = fsp.erode(mg)
mg,_,_ = sp.erode(mg, dt, K_if_used='K_values')
#plot long profiles along channels
if numpy.allclose(elapsed_time%out_tstep,0.) or numpy.allclose(elapsed_time%out_tstep,1.):
pylab.figure("long_profiles")
profile_IDs = prf.channel_nodes(mg, mg.at_node['topographic__steepest_slope'],
mg.at_node['drainage_area'], mg.at_node['flow__receiver_node'])
dists_upstr = prf.get_distances_upstream(mg, len(mg.at_node['topographic__steepest_slope']),
profile_IDs, mg.at_node['flow__link_to_receiver_node'])
prf.plot_profiles(dists_upstr, profile_IDs, mg.at_node['topographic_elevation'])
x_profiles.append(dists_upstr[:])
z_profiles.append(mg.at_node['topographic_elevation'][profile_IDs])
S_profiles.append(mg.at_node['steepest_slope'][profile_IDs])
A_profiles.append(mg.at_node['drainage_area'][profile_IDs])
prf.plot_profiles(dists_upstr, profile_IDs, mg.at_node['topographic__elevation'])
#add uplift
mg.at_node['topographic__elevation'][mg.core_nodes] += 5.*uplift*dt
elapsed_time += dt
precip_perturb = PrecipitationDistribution(input_file=input_file_string, mean_storm=10., mean_interstorm=40., total_t=time_to_run)
out_interval = 5000.
last_trunc = time_to_run #we use this to trigger taking an output plot
for (interval_duration, rainfall_rate) in precip_perturb.yield_storm_interstorm_duration_intensity():
if rainfall_rate != 0.:
mg.at_node['water__unit_flux_in'].fill(rainfall_rate)
mg = fr.run_one_step() #the runoff_rate should pick up automatically
#print 'Area: ', numpy.max(mg.at_node['drainage_area'])
mg,_,_ = sp.erode(mg, interval_duration, Q_if_used='surface_water__discharge', K_if_used='K_values')
#plot long profiles along channels
this_trunc = precip_perturb.elapsed_time//out_interval
if this_trunc != last_trunc: #a new loop round
print('saving a plot at perturbed loop ', out_interval*this_trunc)
pylab.figure("long_profiles")
profile_IDs = prf.channel_nodes(mg, mg.at_node['topographic__steepest_slope'],
mg.at_node['drainage_area'], mg.at_node['flow__receiver_node'])
dists_upstr = prf.get_distances_upstream(mg, len(mg.at_node['topographic__steepest_slope']),
profile_IDs, mg.at_node['flow__link_to_receiver_node'])
prf.plot_profiles(dists_upstr, profile_IDs, mg.at_node['topographic__elevation'])
last_trunc=this_trunc
x_profiles.append(dists_upstr[:])
z_profiles.append(mg.at_node['topographic_elevation'][profile_IDs])
S_profiles.append(mg.at_node['steepest_slope'][profile_IDs])
A_profiles.append(mg.at_node['drainage_area'][profile_IDs])
#add uplift
mg.at_node['topographic__elevation'][mg.core_nodes] += 5.*uplift*interval_duration
#Finalize and plot
elev = mg['node']['topographic__elevation']
elev_r = mg.node_vector_to_raster(elev)
for i in range(nt):
# print 'loop ', i
mg.at_node["topographic__elevation"][mg.core_nodes] += uplift_per_step
mg = fr.route_flow(grid=mg)
if DL_or_TL == "TL":
mg, _ = tle.erode(mg, dt)
else:
mg, _, _ = spe.erode(mg, dt=dt)
if i % init_interval == 0:
print("loop {0}".format(i))
print(
"max_slope",
np.amax(mg.at_node["topographic__steepest_slope"][mg.core_nodes]),
)
pylab.figure("long_profiles_init")
profile_IDs = prf.channel_nodes(
mg,
mg.at_node["topographic__steepest_slope"],
mg.at_node["drainage_area"],
mg.at_node["flow__receiver_node"],
)
dists_upstr = prf.get_distances_upstream(
mg,
len(mg.at_node["topographic__steepest_slope"]),
profile_IDs,
mg.at_node["flow__link_to_receiver_node"],
)
prf.plot_profiles(
dists_upstr, profile_IDs, mg.at_node["topographic__elevation"]
)
print("completed run to steady state...")
else:
mg, _, _ = spe.erode(mg, dt=dt)
if i % init_interval == 0:
print("loop {0}".format(i))
print(
"max_slope",
np.amax(mg.at_node["topographic__steepest_slope"][mg.core_nodes]),
)
pylab.figure("long_profiles_init")
profile_IDs = prf.channel_nodes(
mg,
mg.at_node["topographic__steepest_slope"],
mg.at_node["drainage_area"],
mg.at_node["flow__receiver_node"],
)
dists_upstr = prf.get_distances_upstream(
mg,
len(mg.at_node["topographic__steepest_slope"]),
profile_IDs,
mg.at_node["flow__link_to_receiver_node"],
)
prf.plot_profiles(
dists_upstr, profile_IDs, mg.at_node["topographic__elevation"]
)
print("completed run to steady state...")
if show_figs_in_run:
show() # will cause a hang
# save a copy of the init conditions:
mg_init = deepcopy(mg)
fr = FlowRouter(mg)
sp = SPEroder(mg, input_file)
diffuse = PerronNLDiffuse(mg, input_file)
lin_diffuse = DiffusionComponent(grid=mg, input_stream=input_file)
#perform the loops:
for i in xrange(nt):
#note the input arguments here are not totally standardized between modules
#mg = diffuse.diffuse(mg, i*dt)
mg = lin_diffuse.diffuse(mg, dt)
mg = fr.route_flow(grid=mg)
sp.erode(dt)
##plot long profiles along channels
pylab.figure(6)
profile_IDs = prf.channel_nodes(mg, mg.at_node['steepest_slope'],
mg.at_node['drainage_area'], mg.at_node['upstream_ID_order'],
mg.at_node['flow_receiver'])
dists_upstr = prf.get_distances_upstream(mg, len(mg.at_node['steepest_slope']),
profile_IDs, mg.at_node['links_to_flow_receiver'])
prf.plot_profiles(dists_upstr, profile_IDs, mg.at_node['topographic_elevation'])
print 'Completed loop ', i
print 'Completed the simulation. Plotting...'
#Finalize and plot
# Clear previous plots
pylab.figure(1)
pylab.close()
pylab.figure(1)
sde = SedDepEroder(mg, input_file)
# don't allow overwriting of these, just in case
try:
x_profiles
except NameError:
x_profiles = []
z_profiles = []
S_profiles = []
A_profiles = []
# plot init conds
if make_output_plots:
mg = fr.route_flow(grid=mg)
pylab.figure("long_profile_anim")
ylim([0, y_max])
prf.analyze_channel_network_and_plot(mg)
savefig("0profile_anim_init.png")
close("long_profile_anim")
(profile_IDs, dists_upstr) = prf.analyze_channel_network_and_plot(mg)
start_node = [profile_IDs[0]]
time_on = time()
# perform the loops:
for i in range(nt):
# print 'loop ', i
mg.at_node["topographic__elevation"][mg.core_nodes] += uplift_per_step
mg = fr.run_one_step()
# mg.calc_grad_across_cell_faces(mg.at_node['topographic__elevation'])
# neighbor_slopes = mg.calc_grad_along_node_links(mg.at_node['topographic__elevation'])
# mean_slope = np.mean(np.fabs(neighbor_slopes),axis=1)
# max_slope = np.max(np.fabs(neighbor_slopes),axis=1)