Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def test_fixed_gradient_and_value_boundary():
"""testing multiple boundary conditions with fixed links"""
grid = RasterModelGrid((3, 4))
grid["node"]["topographic__elevation"] = np.zeros(grid.number_of_nodes)
grid["link"]["topographic__slope"] = np.zeros(grid.number_of_links)
grid.set_closed_boundaries_at_grid_edges(False, True, False, True)
grid.status_at_node[4] = FG
grid.status_at_node[7] = FG
grid.status_at_node[1] = 1
assert_array_equal(
grid.status_at_node, [CB, FV, CB, CB, FG, 0, 0, FG, CB, CB, CB, CB]
)
assert_array_equal(
grid.status_at_link, [4, 4, 4, 4, 0, 4, 4, 2, 0, 2, 4, 4, 4, 4, 4, 4, 4]
)
def test_set_status_with_array():
"""Test that active links are reset after changing the node status."""
grid = RasterModelGrid((4, 5))
assert_array_equal(
grid.active_links,
[5, 6, 7, 9, 10, 11, 12, 14, 15, 16, 18, 19, 20, 21, 23, 24, 25],
)
grid.status_at_node[:5] = CB
assert_array_equal(
grid.status_at_node,
[CB, CB, CB, CB, CB, FV, 0, 0, 0, FV, FV, 0, 0, 0, FV, FV, FV, FV, FV, FV],
)
assert_array_equal(
grid.active_links, [9, 10, 11, 12, 14, 15, 16, 18, 19, 20, 21, 23, 24, 25]
)
def test_multiple_status_node(four_by_four_raster):
four_by_four_raster.set_closed_boundaries_at_grid_edges(True, True, False, False)
vals = constant(
four_by_four_raster,
"values",
"node",
where=[CORE_NODE, CLOSED_BOUNDARY],
value=10.0,
)
true_array = np.array(
[
0.0,
0.0,
0.0,
0.0,
0.0,
10.0,
10.0,
10.0,
0.0,
10.0,
10.0,
10.0,
from six.moves import range
import numpy as np
from matplotlib.pyplot import show, plot, figure
from landlab import RasterModelGrid, CLOSED_BOUNDARY, imshow_grid_at_node
from landlab.components import FlowRouter, SedDepEroder
mg = RasterModelGrid((30, 3), 200.) # ((10, 3), 200.)
for edge in (mg.nodes_at_left_edge, mg.nodes_at_top_edge,
mg.nodes_at_right_edge):
mg.status_at_node[edge] = CLOSED_BOUNDARY
z = mg.add_zeros('node', 'topographic__elevation')
th = mg.add_zeros('node', 'channel_sediment__depth')
th += 0.0007
fr = FlowRouter(mg)
sde = SedDepEroder(mg, K_sp=1.e-5,
sed_dependency_type='almost_parabolic',
Qc='power_law', K_t=1.e-5)
z[:] = mg.node_y/10000.
initz = z.copy()
dt = 100.
up = 0.0005
dem_grid.set_nodata_nodes_to_inactive(z, 0) # set nodata nodes to inactive bounds
outlet_node = dem_grid.grid_coords_to_node_id(outlet_row, outlet_column)
grid = RasterModelGrid(4, 5, 1.0)
grid.set_inactive_boundaries(False, True, True, True)
z = grid.add_zeros("node", "Land_Surface__Elevation")
z[6] = 4.5
z[7] = 3.
z[8] = 1.
z[11] = 4.
z[12] = 2.8
z[13] = 2.
# Get array of interior (active) node IDs
interior_nodes = np.where(grid.status_at_node != CLOSED_BOUNDARY)[0]
# Route flow
flow_router = FlowAccumulator(grid, flow_director="D8")
grid = flow_router.route_flow()
for i in range(grid.number_of_nodes):
print(
i,
grid.node_x[i],
grid.node_y[i],
z[i],
grid.status_at_node[i],
r[i],
a[i],
q[i],
ss[i],
zoutlet=z[outlet_loc]
#some helper and parameter values
total_run_time = 500000 #yr
uplift_rate = 0.001 #m/yr
rain_rate = 1 #m/yr
storm_duration = 50 #years
elapsed_time = 0 #years
#set up fluvial incision component
incisor = PowerLawIncision('input_file.txt',rg)
#print "elevations before ", rg.node_vector_to_raster(z)
#interior_nodes are the nodes on which you will be operating
interior_nodes = np.where(rg.status_at_node != CLOSED_BOUNDARY)[0]
while elapsed_time < total_run_time:
#uplift the landscape
z[interior_nodes] = z[interior_nodes]+uplift_rate * storm_duration
z[outlet_loc]=zoutlet
#erode the landscape
z = incisor.run_one_storm(rg,z,rain_rate,storm_duration)
#update the time
elapsed_time = elapsed_time+storm_duration
#below purely for plotting reasons
#instantiate variable of type RouteFlowD8 Class
flow_router = RouteFlowD8(len(z))
#initial flow direction
flowdirs, max_slopes = flow_router.calc_flowdirs(rg,z)
#insantiate variable of type AccumFlow Class
#some helper and parameter values
total_run_time = 10000 #yr
one_twentieth_time = total_run_time/20
uplift_rate = 0.001 #m/yr
rain_rate = 1 #m/yr
storm_duration = 50 #years
elapsed_time = 0 #years
#set up fluvial incision component
incisor = PowerLawIncision('input_file.txt',rg)
#print "elevations before ", rg.node_vector_to_raster(z)
#interior_nodes are the nodes on which you will be operating
interior_nodes = np.where(rg.status_at_node != CLOSED_BOUNDARY)[0]
while elapsed_time < total_run_time:
#erode the landscape
z = incisor.run_one_storm(rg,z,rain_rate,storm_duration)
#uplift the landscape
z[interior_nodes] = z[interior_nodes]+uplift_rate * storm_duration
#update the time
elapsed_time = elapsed_time+storm_duration
#print "elapsed_time", elapsed_time
if elapsed_time%one_twentieth_time == 0:
print("elapsed time",elapsed_time)
elev_raster = rg.node_vector_to_raster(z,True)
# Plot topography
pylab.figure(22)
im = pylab.imshow(elev_raster, cmap=pylab.cm.RdBu, extent=[0, nc*rg.dx, 0, nr*rg.dx])
cb = pylab.colorbar(im)
shape = (mg.number_of_nodes, 1)
plt.quiver(mg.x_of_node.reshape(shape), mg.y_of_node.reshape(shape),
xdist.reshape(shape), ydist.reshape(shape),
color=colors,
angles='xy',
scale_units='xy',
scale=1,
zorder=3)
# Plot differen types of nodes:
o, = plt.plot(mg.x_of_node[mg.status_at_node == CORE_NODE], mg.y_of_node[mg.status_at_node == CORE_NODE], 'b.', label='Core Nodes', zorder=4)
fg, = plt.plot(mg.x_of_node[mg.status_at_node == FIXED_VALUE_BOUNDARY], mg.y_of_node[mg.status_at_node == FIXED_VALUE_BOUNDARY], 'c.', label='Fixed Gradient Nodes', zorder=5)
fv, = plt.plot(mg.x_of_node[mg.status_at_node == FIXED_GRADIENT_BOUNDARY], mg.y_of_node[mg.status_at_node ==FIXED_GRADIENT_BOUNDARY], 'g.', label='Fixed Value Nodes', zorder=6)
c, = plt.plot(mg.x_of_node[mg.status_at_node == CLOSED_BOUNDARY], mg.y_of_node[mg.status_at_node ==CLOSED_BOUNDARY], 'r.', label='Closed Nodes', zorder=7)
node_id = np.arange(mg.number_of_nodes)
flow_to_self = receivers[:,0]==node_id
fts, = plt.plot(mg.x_of_node[flow_to_self], mg.y_of_node[flow_to_self], 'kx', markersize=6, label = 'Flows To Self', zorder=8)
ax = plt.gca()
ax.legend(labels = ['Core Nodes', 'Fixed Gradient Nodes', 'Fixed Value Nodes', 'Closed Nodes', 'Flows To Self'],
handles = [o, fg, fv, c, fts], numpoints=1, loc='center left', bbox_to_anchor=(1.7, 0.5))
sm = plt.cm.ScalarMappable(cmap=propColor, norm=plt.Normalize(vmin=0, vmax=1))
sm._A = []
cx = plt.colorbar(sm)
cx.set_label('Proportion of Flow')
plt.title(title)
h = self.h
g = self.g
n = self.n
tau_crit = self.tau_crit
q = self.q
qs = self.qs
tau = self.tau
rhog = self.rhog
alpha = self.alpha
mpm = self.mpm
zm = self.zm
dtmax = self.dtmax
erode_start_time = self.erode_start_time
ten_thirds = 10./3.
interior_cells = np.where(grid.status_at_node != CLOSED_BOUNDARY)[0]
# Calculate the effective flow depth at active links. Bates et al. 2010
# recommend using the difference between the highest water-surface
# and the highest bed elevation between each pair of cells.
zmax = grid.max_of_link_end_node_values(z)
w = h+z # water-surface height
wmax = grid.max_of_link_end_node_values(w)
hflow = wmax - zmax
# Calculate water-surface slopes
water_surface_slope = grid.calc_grad_of_active_link(w)
# Calculate the unit discharges (Bates et al., eq 11)
q = (q-g*hflow*dtmax*water_surface_slope)/ \
(1.+g*hflow*dtmax*n*n*abs(q)/(hflow**ten_thirds))
elevs = return_array_at_node(grid, elevs)
### Step 1, some basic set-up, gathering information about the grid.
# Calculate the number of nodes.
num_nodes = len(elevs)
# Set the number of receivers and facets.
num_receivers = 2
num_facets = 8
# Create a node array
node_id = np.arange(num_nodes)
# find where there are closed nodes.
closed_nodes = grid.status_at_node == CLOSED_BOUNDARY
# create an array of the triangle numbers
tri_numbers = np.arange(num_facets)
### Step 3, create some triangle datastructures because landlab (smartly)
# makes it hard to deal with diagonals.
# create list of triangle neighbors at node. Use orientation associated
# with tarboton's 1997 algorithm, orthogonal link first, then diagonal.
# has shape, (nnodes, 8 triangles, 2 neighbors)
n_at_node = grid.adjacent_nodes_at_node
dn_at_node = grid.diagonal_adjacent_nodes_at_node
triangle_neighbors_at_node = np.stack(
[
np.vstack((n_at_node[:, 0], dn_at_node[:, 0])),