LinearDiffuser: Model soil creep using “linear diffusion” transport law (no depth dependence)
This section contains documentation and API reference information for the following categories of components:
Two objects based on the EventLayers object exist to make it easier to deal with spatially variable lithology and associated properties. The Lithology components contain information about spatially variable lithology and connect with the Landlab model grid so that when rock is eroded or advected upward by rock uplift the values of rock propeties at the topographic surface are updated.
First is the Lithology component which is a generic object for variable lithology.
Second is LithoLayers which makes it easy to make layered rock.
ChannelProfiler(grid, channel_definition_field='drainage_area', number_of_watersheds=1, minimum_outlet_threshold=0, main_channel_only=True, outlet_nodes=None, minimum_channel_threshold=0, cmap='viridis')[source]¶Bases: landlab.components.profiler.base_profiler._BaseProfiler
Extract and plot the channel profiles in drainage networks.
The ChannelProfiler extracts channel networks from a landlab grid.
In order to extract channel networks, the flow connectivity across the grid must already be identified. This is typically done with the FlowAccumulator component. However, this component does not require that the FlowAccumulator was used. Instead it expects that the following at-node grid fields will be present:
'flow__receiver_node'
'flow__link_to_receiver_node'
The ChannelProfiler can work on grids that have used route-to-one or route-to-multiple flow directing.
To understand how this component works it is useful to define the following terms: watershed, outlet, headwater node, segment.
A watershed is all model grid nodes that drain to a single node, called
the outlet. Channels nodes are identified as nodes that have a
channel_definition_field value greater than or equal to the
minimum_channel_threshold. This channel_definition_field is often
the drainage area (this component’s default). We use a flexible field
rather than only drainage area to support alternative bases for channel
extraction.
The default behaviour of this component is to use an exclusive definition
of a watershed. That is, the two largest watersheds are defined as the
watersheds upstream of the two nodes on the model grid boundary with the
largest values of the channel_definition_field rather than potentially
nested watersheds. Nested watersheds are supported through the use of the
outlet_nodes keyword argument.
Consider the following grid with 10 columns and 7 rows. In this grid is one
watershed with an outlet node indicated by o. Here X indicates
nodes that are not part of the channel network (based on the
channel_definition_field) and . indicates nodes that are part of
the network.
In this and the following examples, we will use only D4 connectivity. The ChannelProfiler, however, knows nothing of connectivity other than what is implied by the two required grid fields.
X X X X X X X X X X
X . X X X X X X X X
X . . X X X . . . X
X X . . X X . X X X
X X X . . . . X X X
X X X . X X X X X X
X X X o X X X X X X
This component can extract the channel network from one or more watersheds.
This option is specified with the keyword argument
number_of_watersheds.
The headwater nodes, shown as @ are nodes that have no upstream
nodes with sufficient area to be classified as a channel.
X X X X X X X X X X
X @ X X X X X X X X
X . . X X X . . @ X
X X . . X X . X X X
X X X . . . . X X X
X X X . X X X X X X
X X X o X X X X X X
For each watershed, the ChannelProfiler either extracts the largest channel
(again, based on the channel_definition_field) or all channel segments
with sufficent values in the channel_definition_field.
Default behavior of this component is to extract only the largest channel in the single largest watershed. This would extract the following channel segment (indicated by the * s).
X X X X X X X X X X
X . X X X X X X X X
X . . X X X * * * X
X X . . X X * X X X
X X X * * * * X X X
X X X * X X X X X X
X X X * X X X X X X
This component verifies that all watershed outlets have a value in the
channel_definition_field of at least minimum_outlet_threshold
(default is 0 units). If no watersheds exist that meet this criteria, an
error is raised.
If a user knows exactly which node or nodes they want to use as the outlet
nodes, then this can be specified using the outlet_nodes keyword
argument. Otherwise the number_of_watersheds (default 1) nodes with the
largest value in the channel_definition_field will be selected as the
outlet nodes from the model grid boundary nodes. Setting
number_of_watersheds to None results in selecting all nodes at the
model grid boundary that meet the criteria for an outlet based on the
channel_definition_field and the minimum_outlet_threshold.
The node IDs and distances upstream of the channel network are stored in
data_structure. It is a dictionary with keys indicating the outlet
node.
For each watershed outlet, the value in the data_structure is itself
a dictionary with keys that are a segment ID tuple of the
(dowstream, upstream) nodes IDs of each channel segment.
For our simple example, these are the node IDs:
X X X X X X X X X X
X 51 X X X X X X X X
X 41 42 X X X 46 47 48 X
X X 32 33 X X 36 X X X
X X X 23 24 25 26 X X X
X X X 13 X X X X X X
X X X 3 X X X X X X
So for our main channel only example, the outlet has an ID of 3, the downstream end of the channel segment is 3, and the upstream end is 48.
The value associated with the segment ID tuple (3, 48) is itself a
dictionary. It has three key-value pairs. First, "ids" contains a list
of the segment node ids ordered from downstream to upstream. It includes
the endpoints. Second, "distances" contains a list of distances
upstream that mirrors the list in "ids". Finally, "color" is an
RGBA tuple indicating the color for the segment.
By default a unique color will be assigned to each watershed. To change the
color, a user can change values stored in data_structure.
Additionally, a cmap keyword argument can provide some user control
over the color at the instantiation of the component.
In the main channel only example, the data structure will look as follows:
{3: {
(3, 48) : {
"ids": [3, 13, 23, 24, 25, 26, 36, 46, 47, 48],
"distances": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
"color": (1, 0, 1, 1),
}
}
}
Three channel segments are idendified if main_channel_only=False.
X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
X . X X X X X X X X X . X X X X X X X X X * X X X X X X X X
X . . X X X . . . X X . . X X X * * * X X * * X X X . . . X
X X . . X X . X X X X X . . X X * X X X X X * * X X . X X X
X X X * . . . X X X X X X * * * * X X X X X X * . . . X X X
X X X * X X X X X X X X X . X X X X X X X X X . X X X X X X
X X X * X X X X X X X X X . X X X X X X X X X . X X X X X X
The data structure associated with this set of three segments is
{3: { (3, 23) : { "ids": [3, 13, 23], "distances": [0, 1, 2] "color": (1, 0, 1, 1), }, (23, 48) : { "ids": [23, 24, 25, 26, 36, 46, 47, 48], "distances": [2, 3, 4, 5, 6, 7, 8, 9], "color": (1, 0, 1, 1), }, (23, 51) : { "ids": [23, 33, 32, 42, 41, 51], "distances": [2, 3, 4, 5, 6, 7] "color": (1, 0, 1, 1), }, } }
Note that the distances upstream are relative to the outlet, not the downstream end of the stream segment.
Next consider a model grid with two watersheds.
X X X X X X X X X X
X . . X X X . X X X
X . X X . X . X X X
o . . . . X . X X X
X X X X X X . X X X
X X X . . . . X X X
X X X X X X . . . X
X X X X X X X X o X
And the following node IDs.
X X X X X X X X X X
X 61 62 X X X 66 X X X
X 51 X X 54 X 56 X X X
40 41 42 43 44 X 46 X X X
X X X X X X 36 X X X
X X X 23 24 25 26 X X X
X X X X X X 16 17 18 X
X X X X X X X X 8 X
The data structure for number_of_watersheds=2 and
main_channel_only=False will be as follows. Note that each watershed
has been assigned a different color tuple value. Here the default viridis
cmap is used.
{8: {
(8, 26) : {
"ids": [8, 18, 17, 16, 26],
"distances": [0, 1, 2, 3, 4]
"color": [ 0.13, 0.57, 0.55, 1. ],
},
(26, 23) : {
"ids": [26, 25, 24, 23],
"distances": [4, 5, 6, 7],
"color": [ 0.13, 0.57, 0.55, 1. ],
},
(26, 66) : {
"ids": [26, 36, 46, 56, 66],
"distances": [4, 5, 6, 7, 8]
"color": [ 0.13, 0.57, 0.55, 1. ],
},
},
40: {
(40, 41) : {
"ids": [40, 41],
"distances": [0, 1]
"color": [ 0.27, 0. , 0.33, 1. ],
},
(41, 54) : {
"ids": [41, 42, 43, 44, 54],
"distances": [2, 3, 4, 5, 6],
"color": [ 0.27, 0. , 0.33, 1. ],
},
(41, 62) : {
"ids": [41, 51, 61, 62],
"distances": [1, 2, 3, 4]
"color": [ 0.27, 0. , 0.33, 1. ],
},
}
}
Examples
Start by importing necessary modules
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator, ChannelProfiler
Create the second example grid we showed above. Note that in order to do this we need to enter the elevations starting from the lower left so the elevation order may seem upside-down. In addition, in this example, elevation is only provided along the profiles. The third line of code below sets all nodes with a value of zero to closed, such that these nodes are igored.
>>> z = np.array([ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
... 0, 0, 0, 0, 0, 0, 4, 3, 2, 0,
... 0, 0, 0, 8, 7, 6, 5, 0, 0, 0,
... 0, 0, 0, 0, 0, 0, 6, 0, 0, 0,
... 1, 3, 4, 5, 6, 0, 7, 0, 0, 0,
... 0, 4, 0, 0, 7, 0, 8, 0, 0, 0,
... 0, 5, 6, 0, 0, 0, 9, 0, 0, 0,
... 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,], dtype=float)
>>> mg = RasterModelGrid((8, 10))
>>> z = mg.add_field("topographic__elevation", z, at="node")
>>> mg.set_nodata_nodes_to_closed(z, 0)
>>> fa = FlowAccumulator(mg, flow_director='D4')
>>> fa.run_one_step()
>>> fa.node_drainage_area.reshape(mg.shape)
array([[ 0., 0., 0., 0., 0., 0., 0., 0., 11., 0.],
[ 0., 0., 0., 0., 0., 0., 9., 10., 11., 0.],
[ 0., 0., 0., 1., 2., 3., 8., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 4., 0., 0., 0.],
[ 8., 8., 4., 3., 2., 0., 3., 0., 0., 0.],
[ 0., 3., 0., 0., 1., 0., 2., 0., 0., 0.],
[ 0., 2., 1., 0., 0., 0., 1., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])
>>> profiler = ChannelProfiler(
... mg,
... number_of_watersheds=2,
... minimum_channel_threshold=0,
... main_channel_only=False)
>>> profiler.run_one_step()
The keys of the property data_structure are the IDs of the two
outlet nodes.
>>> profiler.data_structure.keys()
odict_keys([40, 8])
Within the data structure, the value at key 40, is a dictionary of the
three segments, each specified by a (dowstream, upstream) tuple:
>>> profiler.data_structure[40].keys()
odict_keys([(40, 41), (41, 54), (41, 62)])
The value of the segment between nodes 40 and 41 has the following components:
>>> profiler.data_structure[40][(40, 41)]["ids"]
array([40, 41])
>>> profiler.data_structure[40][(40, 41)]["distances"]
array([ 0., 1.])
>>> np.round(profiler.data_structure[40][(40, 41)]["color"], decimals=2)
array([ 0.27, 0. , 0.33, 1. ])
A parallel structure exists for the segment between nodes 41 and 54:
>>> profiler.data_structure[40][(41, 54)]["ids"]
array([41, 42, 43, 44, 54])
>>> profiler.data_structure[40][(41, 54)]["distances"]
array([ 1., 2., 3., 4., 5.])
>>> np.round(profiler.data_structure[40][(41, 54)]["color"], decimals=2)
array([ 0.27, 0. , 0.33, 1. ])
And the segment between nodes 41 and 62.
>>> profiler.data_structure[40][(41, 62)]["ids"]
array([41, 51, 61, 62])
>>> profiler.data_structure[40][(41, 62)]["distances"]
array([ 1., 2., 3., 4.])
>>> np.round(profiler.data_structure[40][(41, 62)]["color"], decimals=2)
array([ 0.27, 0. , 0.33, 1. ])
The rest of the profile_structure encodes information about the second
watershed, which drains to node 8.
>>> profiler.data_structure[8].keys()
odict_keys([(8, 26), (26, 23), (26, 66)])
>>> profiler.data_structure[8][(8, 26)]["ids"]
array([ 8, 18, 17, 16, 26])
>>> profiler.data_structure[8][(8, 26)]["distances"]
array([ 0., 1., 2., 3., 4.])
>>> np.round(profiler.data_structure[8][(8, 26)]["color"], decimals=2)
array([ 0.13, 0.57, 0.55, 1. ])
>>> profiler.data_structure[8][(26, 23)]["ids"]
array([26, 25, 24, 23])
>>> profiler.data_structure[8][(26, 23)]["distances"]
array([ 4., 5., 6., 7.])
>>> np.round(profiler.data_structure[8][(26, 23)]["color"], decimals=2)
array([ 0.13, 0.57, 0.55, 1. ])
>>> profiler.data_structure[8][(26, 66)]["ids"]
array([26, 36, 46, 56, 66])
>>> profiler.data_structure[8][(26, 66)]["distances"]
array([ 4., 5., 6., 7., 8.])
>>> np.round(profiler.data_structure[8][(26, 66)]["color"], decimals=2)
array([ 0.13, 0.57, 0.55, 1. ])
The ChannelProfiler is designed to be flexible, and by careful combination
of its instantiation variables can be used to extract many useful forms of
profile. In these examples, we will use the default
channel_definition_field, the drainage area.
To illustrate, lets start by creating a landscape model.
>>> from landlab.components import FastscapeEroder
>>> mg = RasterModelGrid((100, 120), xy_spacing=2)
>>> np.random.seed(42)
>>> z = mg.add_zeros('topographic__elevation', at='node')
>>> z[mg.core_nodes] += np.random.randn(mg.core_nodes.size)
>>> fa = FlowAccumulator(mg)
>>> sp = FastscapeEroder(mg, K_sp=0.0001)
>>> dt = 1000
>>> for i in range(200):
... fa.run_one_step()
... sp.run_one_step(dt=dt)
... z[mg.core_nodes] += 0.001 * dt
Some options:
Default: Extract a the single biggest channel draining to the model grid boundary traced back all the way to the watershed divide.
>>> profiler = ChannelProfiler(mg)
Extract the largest channel draining to each of the four largest outlet nodes on the model grid boundary traced back all the way to the watershed divide.
>>> profiler = ChannelProfiler(mg, number_of_watersheds=4)
Extract the single largest channel draining to node 2933. Note that the
keyword argument outlet_nodes must be an iterable.
>>> profiler = ChannelProfiler(mg, outlet_nodes=[2933])
Extract the largest channel draining to each of the four largest outlet
nodes on the model grid boundary traced back to nodes with
channel_definition_field values of 500.
>>> profiler = ChannelProfiler(
... mg,
... number_of_watersheds=4,
... minimum_channel_threshold=500)
Extract a the single biggest channel draining to the model grid boundary
based on the field surface_water__discharge traced back to discharge
values of 500.
>>> profiler = ChannelProfiler(mg,
... channel_definition_field='surface_water__discharge',
... minimum_channel_threshold=500)
Extract the single largest channel within all watersheds with an outlet
with channel_definition_field greater than 1e3. Trace the channels
up to the point in each watershed in which the channels have values in the
channel_definition_field of 500.
>>> profiler = ChannelProfiler(
... mg,
... number_of_watersheds=None,
... minimum_outlet_threshold=1e3,
... minimum_channel_threshold=500)
Extract two trunk channels beginning at the given nodes, traced up to a
a minimum channel_definition_field value of of 500. Note that
number_of_watersheds must match the size of outlet_nodes.
>>> profiler = ChannelProfiler(
... mg,
... outlet_nodes=[6661, 6250],
... number_of_watersheds=2,
... minimum_channel_threshold=500)
Extract every possible channel (not just the largest one), leading from the
four highest model grid boundary nodes traced back to a
channel_definition_field threshold of 20.
>>> profiler = ChannelProfiler(mg,
... number_of_watersheds=4,
... main_channel_only=False,
... minimum_channel_threshold=20)
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
None Listed
| Parameters: |
|
|---|
__init__(grid, channel_definition_field='drainage_area', number_of_watersheds=1, minimum_outlet_threshold=0, main_channel_only=True, outlet_nodes=None, minimum_channel_threshold=0, cmap='viridis')[source]¶| Parameters: |
|
|---|
assign_colors(color_mapping=None)[source]¶Assign a unique color for each watershed.
| Parameters: | color_mapping (str) – Color map name. |
|---|
data_structure¶OrderedDict defining the channel network.
The IDs and upstream distance of the channel network nodes are stored
in data_structure. It is a dictionary with keys of the outlet node
ID.
For each watershed outlet, the value in the data_structure is
itself a dictionary with keys that are a segment ID tuple of the
(dowstream, upstream) nodes IDs of each channel segment.
The value associated with the segment ID tuple
(dowstream, upstream) is itself a dictionary. It has three
key-value pairs. First, "ids" contains a list of the segment node
IDs ordered from downstream to upstream. It includes the endpoints.
Second, "distances" contains a list of distances upstream that
mirrors the list in "ids". Finally, "color" is an RGBA tuple
indicating the color for the segment.
ChiFinder(grid, reference_concavity=0.5, min_drainage_area=1000000.0, reference_area=1.0, use_true_dx=False, clobber=False)[source]¶Bases: landlab.core.model_component.Component
Calculate Chi Indices.
This component calculates chi indices, sensu Perron & Royden, 2013, for a Landlab landscape.
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator, FastscapeEroder
>>> from landlab.components import ChiFinder
>>> mg = RasterModelGrid((3, 4))
>>> for nodes in (mg.nodes_at_right_edge, mg.nodes_at_bottom_edge,
... mg.nodes_at_top_edge):
... mg.status_at_node[nodes] = mg.BC_NODE_IS_CLOSED
>>> _ = mg.add_field("topographic__elevation", mg.node_x, at="node")
>>> fr = FlowAccumulator(mg, flow_director='D8')
>>> cf = ChiFinder(mg,
... min_drainage_area=1.,
... reference_concavity=1.)
>>> fr.run_one_step()
>>> cf.calculate_chi()
>>> mg.at_node['channel__chi_index'].reshape(mg.shape)[1, :]
array([ 0.5, 1. , 2. , 0. ])
>>> mg2 = RasterModelGrid((5, 5), xy_spacing=100.)
>>> for nodes in (mg2.nodes_at_right_edge, mg2.nodes_at_bottom_edge,
... mg2.nodes_at_top_edge):
... mg2.status_at_node[nodes] = mg2.BC_NODE_IS_CLOSED
>>> _ = mg2.add_zeros('node', 'topographic__elevation')
>>> mg2.at_node['topographic__elevation'][mg2.core_nodes] = mg2.node_x[
... mg2.core_nodes]/1000.
>>> np.random.seed(0)
>>> mg2.at_node['topographic__elevation'][
... mg2.core_nodes] += np.random.rand(mg2.number_of_core_nodes)
>>> fr2 = FlowAccumulator(mg2, flow_director='D8')
>>> sp2 = FastscapeEroder(mg2, K_sp=0.01)
>>> cf2 = ChiFinder(
... mg2,
... min_drainage_area=0.,
... reference_concavity=0.5)
>>> for i in range(10):
... mg2.at_node['topographic__elevation'][mg2.core_nodes] += 10.
... fr2.run_one_step()
... sp2.run_one_step(1000.)
>>> fr2.run_one_step()
>>> cf2.calculate_chi()
>>> mg2.at_node['channel__chi_index'].reshape(
... mg2.shape) # doctest: +NORMALIZE_WHITESPACE
array([[ 0. , 0. , 0. , 0. , 0. ],
[ 0.77219416, 1.54438833, 2.63643578, 2.61419437, 0. ],
[ 1.09204746, 2.18409492, 1.52214691, 2.61419437, 0. ],
[ 0.44582651, 0.89165302, 1.66384718, 2.75589464, 0. ],
[ 0. , 0. , 0. , 0. , 0. ]])
>>> cf3 = ChiFinder(
... mg2,
... min_drainage_area=20000.,
... use_true_dx=True,
... reference_concavity=0.5,
... reference_area=mg2.at_node['drainage_area'].max(),
... clobber=True)
>>> cf3.calculate_chi()
>>> cf3.chi_indices.reshape(mg2.shape) # doctest: +NORMALIZE_WHITESPACE
array([[ 0. , 0. , 0. , 0. , 0. ],
[ 0. , 173.20508076, 0. , 0. , 0. ],
[ 0. , 0. , 270.71067812, 0. , 0. ],
[ 0. , 100. , 236.60254038, 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ]])
>>> cf3.hillslope_mask.reshape(mg2.shape)
array([[ True, True, True, True, True],
[False, False, True, True, True],
[ True, True, False, True, True],
[False, False, False, True, True],
[ True, True, True, True, True]], dtype=bool)
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
Perron, J., Royden, L. (2012). An integral approach to bedrock river profile analysis Earth Surface Processes and Landforms 38(6), 570-576. https://dx.doi.org/10.1002/esp.3302
| Parameters: |
|
|---|
__init__(grid, reference_concavity=0.5, min_drainage_area=1000000.0, reference_area=1.0, use_true_dx=False, clobber=False)[source]¶| Parameters: |
|
|---|
best_fit_chi_elevation_gradient_and_intercept(ch_nodes=None)[source]¶Returns least squares best fit for a straight line through a chi plot.
| Parameters: | ch_nodes (array of ints or None) – Nodes at which to consider chi and elevation values. If None, will use all nodes in grid with area greater than the component min_drainage_area. |
|---|---|
| Returns: | coeffs – A len-2 array containing the m then z0, where z = z0 + m * chi. |
| Return type: | array(gradient, intercept) |
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> from landlab.components import ChiFinder
>>> mg = RasterModelGrid((3, 4))
>>> for nodes in (mg.nodes_at_right_edge, mg.nodes_at_bottom_edge,
... mg.nodes_at_top_edge):
... mg.status_at_node[nodes] = mg.BC_NODE_IS_CLOSED
>>> z = mg.add_field("topographic__elevation", mg.node_x.copy(), at="node")
>>> z[4:8] = np.array([0.5, 1., 2., 0.])
>>> fr = FlowAccumulator(mg, flow_director='D8')
>>> cf = ChiFinder(
... mg,
... min_drainage_area=1.,
... reference_concavity=1.)
>>> fr.run_one_step()
>>> cf.calculate_chi()
>>> mg.at_node['channel__chi_index'].reshape(mg.shape)[1, :]
array([ 0.5, 1. , 2. , 0. ])
>>> coeffs = cf.best_fit_chi_elevation_gradient_and_intercept()
>>> np.allclose(np.array([1., 0.]), coeffs)
True
calculate_chi()[source]¶Calculate local chi indices.
This is the main method. Call it to calculate local chi indices at all points with drainage areas greater than min_drainage_area.
Chi of any node without a defined value is reported as 0. These nodes
are also identified in the mask retrieved with hillslope_mask.
chi_indices¶Return the array of channel steepness indices.
Nodes not in the channel receive zeros.
create_chi_plot(channel_heads=None, label_axes=True, symbol='kx', plot_line=False, line_symbol='r-')[source]¶Plots a “chi plot” (chi vs elevation for points in channel network).
If channel_heads is provided, only the channel nodes downstream of the provided points (and with area > min_drainage_area) will be plotted.
| Parameters: |
|
|---|
hillslope_mask¶Return a boolean array, False where steepness indices exist.
integrate_chi_avg_dx(valid_upstr_order, chi_integrand, chi_array, mean_dx)[source]¶Calculates chi at each channel node by summing chi_integrand.
This method assumes a uniform, mean spacing between nodes. Method is deliberately split out for potential cythonization at a later stage.
| Parameters: |
|
|---|
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> from landlab.components import ChiFinder
>>> mg = RasterModelGrid((5, 4))
>>> for nodes in (mg.nodes_at_right_edge, mg.nodes_at_bottom_edge,
... mg.nodes_at_top_edge):
... mg.status_at_node[nodes] = mg.BC_NODE_IS_CLOSED
>>> z = mg.node_x.copy()
>>> z[[5, 13]] = z[6] # guard nodes
>>> _ = mg.add_field("topographic__elevation", z, at="node")
>>> fr = FlowAccumulator(mg, flow_director='D8')
>>> cf = ChiFinder(mg)
>>> fr.run_one_step()
>>> ch_nodes = np.array([4, 8, 12, 5, 9, 13, 6, 10, 14])
>>> ch_integrand = 3.*np.ones(9, dtype=float) # to make calc clearer
>>> chi_array = np.zeros(mg.number_of_nodes, dtype=float)
>>> cf.integrate_chi_avg_dx(ch_nodes, ch_integrand, chi_array, 0.5)
>>> chi_array.reshape(mg.shape)
array([[ 0. , 0. , 0. , 0. ],
[ 1.5, 3. , 4.5, 0. ],
[ 1.5, 3. , 4.5, 0. ],
[ 1.5, 3. , 4.5, 0. ],
[ 0. , 0. , 0. , 0. ]])
integrate_chi_each_dx(valid_upstr_order, chi_integrand_at_nodes, chi_array)[source]¶Calculates chi at each channel node by summing chi_integrand*dx.
This method accounts explicitly for spacing between each node. Method is deliberately split out for potential cythonization at a later stage. Uses a trapezium integration method.
| Parameters: |
|
|---|
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> from landlab.components import ChiFinder
>>> mg = RasterModelGrid((5, 4), xy_spacing=3.)
>>> for nodes in (mg.nodes_at_right_edge, mg.nodes_at_bottom_edge,
... mg.nodes_at_top_edge):
... mg.status_at_node[nodes] = mg.BC_NODE_IS_CLOSED
>>> z = mg.node_x.copy()
>>> z[[5, 13]] = z[6] # guard nodes
>>> _ = mg.add_field("topographic__elevation", z, at="node")
>>> fr = FlowAccumulator(mg, flow_director='D8')
>>> cf = ChiFinder(mg)
>>> fr.run_one_step()
>>> ch_nodes = np.array([4, 8, 12, 5, 9, 13, 6, 10, 14])
>>> ch_integrand = 2.*np.ones(mg.number_of_nodes,
... dtype=float) # to make calc clearer
>>> chi_array = np.zeros(mg.number_of_nodes, dtype=float)
>>> cf.integrate_chi_each_dx(ch_nodes, ch_integrand, chi_array)
>>> chi_array.reshape(mg.shape)
array([[ 0. , 0. , 0. , 0. ],
[ 0. , 6. , 14.48528137, 0. ],
[ 0. , 6. , 12. , 0. ],
[ 0. , 6. , 14.48528137, 0. ],
[ 0. , 0. , 0. , 0. ]])
>>> from landlab.components import FastscapeEroder
>>> mg2 = RasterModelGrid((5, 5), xy_spacing=100.)
>>> for nodes in (mg2.nodes_at_right_edge, mg2.nodes_at_bottom_edge,
... mg2.nodes_at_top_edge):
... mg2.status_at_node[nodes] = mg2.BC_NODE_IS_CLOSED
>>> _ = mg2.add_zeros('node', 'topographic__elevation')
>>> mg2.at_node['topographic__elevation'][mg2.core_nodes] = mg2.node_x[
... mg2.core_nodes]/1000.
>>> np.random.seed(0)
>>> mg2.at_node['topographic__elevation'][
... mg2.core_nodes] += np.random.rand(mg2.number_of_core_nodes)
>>> fr2 = FlowAccumulator(mg2, flow_director='D8')
>>> sp2 = FastscapeEroder(mg2, K_sp=0.01)
>>> cf2 = ChiFinder(mg2, min_drainage_area=1., reference_concavity=0.5,
... use_true_dx=True)
>>> for i in range(10):
... mg2.at_node['topographic__elevation'][mg2.core_nodes] += 10.
... fr2.run_one_step()
... sp2.run_one_step(1000.)
>>> fr2.run_one_step()
>>> output_array = np.zeros(25, dtype=float)
>>> cf2.integrate_chi_each_dx(mg2.at_node['flow__upstream_node_order'],
... np.ones(25, dtype=float),
... output_array)
>>> output_array.reshape(mg2.shape)
array([[ 0. , 0. , 0. , 0. , 0. ],
[ 0. , 100. , 200. , 382.84271247, 0. ],
[ 0. , 100. , 241.42135624, 341.42135624, 0. ],
[ 0. , 100. , 200. , 300. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ]])
masked_chi_indices¶Returns a masked array version of the ‘channel__chi_index’ field. This enables easier plotting of the values with.
landlab.imshow_grid_at_node or similar.
Examples
Make a topographic map with an overlay of chi values:
>>> from landlab import imshow_grid_at_node
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator, FastscapeEroder
>>> from landlab.components import ChiFinder
>>> mg = RasterModelGrid((5, 5), xy_spacing=100.)
>>> for nodes in (mg.nodes_at_right_edge, mg.nodes_at_bottom_edge,
... mg.nodes_at_top_edge):
... mg.status_at_node[nodes] = mg.BC_NODE_IS_CLOSED
>>> _ = mg.add_zeros('node', 'topographic__elevation')
>>> mg.at_node['topographic__elevation'][mg.core_nodes] = mg.node_x[
... mg.core_nodes]/1000.
>>> np.random.seed(0)
>>> mg.at_node['topographic__elevation'][
... mg.core_nodes] += np.random.rand(mg.number_of_core_nodes)
>>> fr = FlowAccumulator(mg, flow_director='D8')
>>> sp = FastscapeEroder(mg, K_sp=0.01)
>>> cf = ChiFinder(mg, min_drainage_area=20000.)
>>> for i in range(10):
... mg.at_node['topographic__elevation'][mg.core_nodes] += 10.
... fr.run_one_step()
... sp.run_one_step(1000.)
>>> fr.run_one_step()
>>> cf.calculate_chi()
>>> imshow_grid_at_node(mg, 'topographic__elevation',
... allow_colorbar=False)
>>> imshow_grid_at_node(mg, cf.masked_chi_indices,
... color_for_closed=None, cmap='winter')
mean_channel_node_spacing(ch_nodes)[source]¶Calculates the mean spacing between all adjacent channel nodes.
| Parameters: | ch_nodes (array of ints) – The nodes within the defined channel network. |
|---|---|
| Returns: | mean_spacing – The mean spacing between all nodes in the network. |
| Return type: | float (m) |
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> from landlab.components import ChiFinder
>>> mg = RasterModelGrid((5, 4), xy_spacing=2.)
>>> for nodes in (mg.nodes_at_right_edge, mg.nodes_at_bottom_edge,
... mg.nodes_at_top_edge):
... mg.status_at_node[nodes] = mg.BC_NODE_IS_CLOSED
>>> z = mg.node_x.copy()
>>> z[[5, 13]] = z[6] # guard nodes
>>> _ = mg.add_field("topographic__elevation", z, at="node")
>>> fr = FlowAccumulator(mg, flow_director='D8')
>>> cf = ChiFinder(mg)
>>> fr.run_one_step()
>>> ch_nodes = np.array([4, 8, 12, 5, 9, 13, 6, 10, 14])
>>> cf.mean_channel_node_spacing(ch_nodes)
2.2761423749153966
nodes_downstream_of_channel_head(channel_head)[source]¶Find and return an array with nodes downstream of channel_head.
| Parameters: | channel_head (int) – Node ID of channel head from which to get downstream nodes. |
|---|
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator, ChiFinder
>>> mg = RasterModelGrid((3, 4))
>>> for nodes in (mg.nodes_at_right_edge, mg.nodes_at_bottom_edge,
... mg.nodes_at_top_edge):
... mg.status_at_node[nodes] = mg.BC_NODE_IS_CLOSED
>>> z = mg.add_field("topographic__elevation", mg.node_x.copy(), at="node")
>>> z[4:8] = np.array([0.5, 1., 2., 0.])
>>> fr = FlowAccumulator(mg, flow_director='D8')
>>> fr.run_one_step()
>>> mg.at_node['flow__receiver_node']
array([ 0, 1, 2, 3, 4, 4, 5, 7, 8, 9, 10, 11])
>>> cf = ChiFinder(mg, min_drainage_area=0., reference_concavity=1.)
>>> cf.calculate_chi()
>>> cf.nodes_downstream_of_channel_head(6)
[6, 5, 4]
DepressionFinderAndRouter(grid, routing='D8', pits='flow__sink_flag', reroute_flow=True)[source]¶Bases: landlab.core.model_component.Component
Find depressions on a topographic surface.
This component identifies depressions in a topographic surface, finds an outlet for each depression. If directed to do so (default True), and the component is able to find existing routing fields output from the ‘route_flow_dn’ component, it will then modify the drainage directions and accumulations already stored in the grid to route flow across these depressions.
Note that in general properties of this class named “depression” identify each individual pit in the topography, including those that will merge once the fill is performed. Those named “lake” return the unique lakes created by the fill, and are probably the properties most users will want.
Note also that the structure of drainage within the lakes is not guaranteed, and in particular, may not be symmetrical even if your boundary conditions are. However, the outputs from the lake will all still be correct.
Note the routing part of this component is not yet compatible with irregular grids.
The prinary method of this class is map_depressions().
Examples
Route flow across a depression in a sloped surface.
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator, DepressionFinderAndRouter
>>> mg = RasterModelGrid((7, 7), xy_spacing=0.5)
>>> z = mg.add_field("topographic__elevation", mg.node_x.copy(), at="node")
>>> z += 0.01 * mg.node_y
>>> mg.at_node['topographic__elevation'].reshape(mg.shape)[2:5, 2:5] *= 0.1
>>> fr = FlowAccumulator(mg, flow_director='D8')
>>> fr.run_one_step() # the flow "gets stuck" in the hole
>>> mg.at_node['flow__receiver_node'].reshape(mg.shape)
array([[ 0, 1, 2, 3, 4, 5, 6],
[ 7, 7, 16, 17, 18, 18, 13],
[14, 14, 16, 16, 17, 18, 20],
[21, 21, 16, 23, 24, 25, 27],
[28, 28, 23, 30, 31, 32, 34],
[35, 35, 30, 31, 32, 32, 41],
[42, 43, 44, 45, 46, 47, 48]])
>>> mg.at_node['drainage_area'].reshape(mg.shape)
array([[ 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
[ 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0. ],
[ 0.25, 0.25, 5. , 1.5 , 1. , 0.25, 0. ],
[ 0.25, 0.25, 3. , 0.75, 0.5 , 0.25, 0. ],
[ 0.25, 0.25, 2. , 1.5 , 1. , 0.25, 0. ],
[ 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. , 0. ]])
>>> df = DepressionFinderAndRouter(mg)
>>> df.map_depressions() # reroute_flow defaults to True
>>> mg.at_node['flow__receiver_node'].reshape(mg.shape)
array([[ 0, 1, 2, 3, 4, 5, 6],
[ 7, 7, 16, 17, 18, 18, 13],
[14, 14, 8, 16, 17, 18, 20],
[21, 21, 16, 16, 24, 25, 27],
[28, 28, 23, 24, 24, 32, 34],
[35, 35, 30, 31, 32, 32, 41],
[42, 43, 44, 45, 46, 47, 48]])
>>> mg.at_node['drainage_area'].reshape(mg.shape)
array([[ 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
[ 5.25, 5.25, 0.25, 0.25, 0.25, 0.25, 0. ],
[ 0.25, 0.25, 5. , 1.5 , 1. , 0.25, 0. ],
[ 0.25, 0.25, 0.75, 2.25, 0.5 , 0.25, 0. ],
[ 0.25, 0.25, 0.5 , 0.5 , 1. , 0.25, 0. ],
[ 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. , 0. ]])
>>> df.lake_at_node.reshape(mg.shape) # doctest: +NORMALIZE_WHITESPACE
array([[False, False, False, False, False, False, False],
[False, False, False, False, False, False, False],
[False, False, True, True, True, False, False],
[False, False, True, True, True, False, False],
[False, False, True, True, True, False, False],
[False, False, False, False, False, False, False],
[False, False, False, False, False, False, False]], dtype=bool)
>>> df.lake_map.reshape(mg.shape) # doctest: +NORMALIZE_WHITESPACE
array([[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, 16, 16, 16, -1, -1],
[-1, -1, 16, 16, 16, -1, -1],
[-1, -1, 16, 16, 16, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1]])
>>> df.lake_codes # a unique code for each lake present on the grid
array([16])
>>> df.lake_outlets # the outlet node of each lake in lake_codes
array([8])
>>> df.lake_areas # the area of each lake in lake_codes
array([ 2.25])
Because rereoute_flow defaults to True, the flow connectivity fields
created by the FlowAccumulator
will have now been modified to route flow over the depressions in the
surface. The topogrphy itself is not modified.
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
Tucker, G. E., Lancaster, S. T., Gasparini, N. M., and Bras, R. L.: The Channel-Hillslope Integrated Landscape Development Model (CHILD), in: Landscape Erosion and Evolution Modeling, Springer US, Boston, MA, USA, 349–388, 2001.
Create a DepressionFinderAndRouter.
Constructor assigns a copy of the grid, sets the current time, and calls the initialize method.
| Parameters: |
|
|---|
__init__(grid, routing='D8', pits='flow__sink_flag', reroute_flow=True)[source]¶Create a DepressionFinderAndRouter.
Constructor assigns a copy of the grid, sets the current time, and calls the initialize method.
| Parameters: |
|
|---|
assign_outlet_receiver(outlet_node)[source]¶Find drainage direction for outlet_node that does not flow into its own lake.
Examples
>>> import numpy as np
>>> from landlab.components import DepressionFinderAndRouter
>>> from landlab.components import FlowAccumulator
>>> from landlab import RasterModelGrid
>>> rg = RasterModelGrid((7, 7))
>>> rg.status_at_node[rg.nodes_at_right_edge] = rg.BC_NODE_IS_CLOSED
>>> z = rg.add_zeros("topographic__elevation", at="node")
>>> z[:] = rg.x_of_node + 0.01 * rg.y_of_node
>>> lake_nodes = np.array([10, 16, 17, 18, 24, 32, 33, 38, 40])
>>> z[lake_nodes] *= 0.1
>>> fr = FlowAccumulator(rg, flow_director='D4')
>>> fr.run_one_step()
>>> rg.at_node['flow__receiver_node']
array([ 0, 1, 2, 3, 4, 5, 6, 7, 7, 16, 10, 10, 11, 13, 14, 14, 16,
16, 17, 18, 20, 21, 21, 16, 17, 24, 33, 27, 28, 28, 29, 24, 32, 32,
34, 35, 35, 38, 38, 38, 33, 41, 42, 43, 44, 45, 46, 47, 48])
>>> df = DepressionFinderAndRouter(rg, routing='D4')
>>> df.map_depressions()
>>> rg.at_node['flow__receiver_node']
array([ 0, 1, 2, 3, 4, 5, 6, 7, 7, 16, 17, 10, 11, 13, 14, 14, 15,
16, 17, 18, 20, 21, 21, 16, 17, 24, 33, 27, 28, 28, 29, 38, 31, 32,
34, 35, 35, 36, 37, 38, 33, 41, 42, 43, 44, 45, 46, 47, 48])
>>> fr = FlowAccumulator(rg, flow_director='D8')
>>> fr.run_one_step()
>>> rg.at_node['flow__receiver_node']
array([ 0, 1, 2, 3, 4, 5, 6, 7, 7, 16, 16, 10, 18, 13, 14, 14, 16,
16, 17, 18, 20, 21, 21, 16, 16, 24, 33, 27, 28, 28, 24, 24, 24, 32,
34, 35, 35, 38, 38, 38, 32, 41, 42, 43, 44, 45, 46, 47, 48])
>>> df = DepressionFinderAndRouter(rg, routing='D8')
>>> df.map_depressions()
>>> rg.at_node['flow__receiver_node']
array([ 0, 1, 2, 3, 4, 5, 6, 7, 7, 16, 16, 10, 18, 13, 14, 14, 8,
16, 17, 18, 20, 21, 21, 16, 16, 24, 33, 27, 28, 28, 24, 24, 24, 32,
34, 35, 35, 38, 32, 38, 32, 41, 42, 43, 44, 45, 46, 47, 48])
depression_depth¶At node array of depression depths.
depression_outlet_map¶At node array indicating the node-id of the depression outlet.
find_depression_from_pit(pit_node, reroute_flow=True)[source]¶Find the extent of the nodes that form a pit.
Identify extent of depression/lake whose lowest point is the node pit_node (which is a itself a pit, a.k.a., closed depression).
| Parameters: | pit_node (int) – The node that is the lowest point of a pit. |
|---|
flood_status¶Map of flood status (_PIT, _CURRENT_LAKE, _UNFLOODED, or _FLOODED).
is_pit¶At node array indicating whether the node is a pit or not.
is_valid_outlet(the_node)[source]¶Check if a node is a valid outlet for the depression.
| Parameters: |
|
|---|---|
| Returns: |
|
| Return type: | boolean |
lake_areas¶A nlakes-long array of the area of each lake.
The order is the same as that returned by lake_codes.
lake_at_node¶Return a boolean array, True if the node is flooded, False otherwise.
lake_codes¶Returns the unique code assigned to each unique lake.
These are the values used to map the lakes in the property “lake_map”.
lake_map¶Return an array of ints, where each node within a lake is labelled with a unique (non-consecutive) code corresponding to each unique lake.
The codes used can be obtained with lake_codes. Nodes not in a lake are labelled with self._grid.BAD_INDEX
lake_outlets¶Returns the unique outlets for each lake, in same order as the return from lake_codes.
lake_volumes¶A nlakes-long array of the volume of each lake.
The order is the same as that returned by lake_codes.
map_depressions()[source]¶Map depressions/lakes in a topographic surface.
Examples
Test #1: 5x5 raster grid with a diagonal lake.
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import (
... DepressionFinderAndRouter)
>>> rg = RasterModelGrid((5, 5))
>>> z = rg.add_zeros("topographic__elevation", at="node")
>>> z[:] = np.array([100., 100., 95., 100., 100.,
... 100., 101., 92., 1., 100.,
... 100., 101., 2., 101., 100.,
... 100., 3., 101., 101., 100.,
... 90., 95., 100., 100., 100.])
>>> df = DepressionFinderAndRouter(rg, reroute_flow=False)
>>> df.map_depressions()
>>> df.display_depression_map() # doctest: +NORMALIZE_WHITESPACE
. . . . .
. . . ~ .
. . ~ . .
. ~ . . .
o . . . .
node_can_drain(the_node)[source]¶Check if a node has drainage away from the current lake/depression.
| Parameters: |
|
|---|---|
| Returns: |
|
| Return type: | boolean |
number_of_lakes¶Return the number of individual lakes.
number_of_pits¶The number of pits on the grid.
pit_node_ids¶Node IDs of grid nodes identified as pits.
receivers¶At node array indicating which node receives flow.
DepthDependentDiffuser(grid, linear_diffusivity=1.0, soil_transport_decay_depth=1.0)[source]¶Bases: landlab.core.model_component.Component
This component implements a depth and slope dependent linear diffusion rule in the style of Johnstone and Hilley (2014).
Hillslope sediment flux uses depth dependent component inspired by Johnstone and Hilley (2014). The flux \(q_s\) is given as:
where \(D\) is is the diffusivity, \(S\) is the slope (defined as negative downward), \(H\) is the soil depth on links, and \(H^*\) is the soil transport decay depth.
This component will ignore soil thickness located at non-core nodes.
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import ExponentialWeatherer
>>> from landlab.components import DepthDependentDiffuser
>>> mg = RasterModelGrid((5, 5))
>>> soilTh = mg.add_zeros("soil__depth", at="node")
>>> z = mg.add_zeros("topographic__elevation", at="node")
>>> BRz = mg.add_zeros("bedrock__elevation", at="node")
>>> expweath = ExponentialWeatherer(mg)
>>> DDdiff = DepthDependentDiffuser(mg)
>>> expweath.calc_soil_prod_rate()
>>> np.allclose(mg.at_node['soil_production__rate'][mg.core_nodes], 1.)
True
>>> DDdiff.run_one_step(2.)
>>> np.allclose(mg.at_node['topographic__elevation'][mg.core_nodes], 0.)
True
>>> np.allclose(mg.at_node['bedrock__elevation'][mg.core_nodes], -2.)
True
>>> np.allclose(mg.at_node['soil__depth'][mg.core_nodes], 2.)
True
Now with a slope:
>>> mg = RasterModelGrid((3, 5))
>>> soilTh = mg.add_zeros("soil__depth", at="node")
>>> z = mg.add_zeros("topographic__elevation", at="node")
>>> BRz = mg.add_zeros("bedrock__elevation", at="node")
>>> z += mg.node_x.copy()
>>> BRz += mg.node_x/2.
>>> soilTh[:] = z - BRz
>>> expweath = ExponentialWeatherer(mg)
>>> DDdiff = DepthDependentDiffuser(mg)
>>> expweath.calc_soil_prod_rate()
>>> np.allclose(
... mg.at_node['soil_production__rate'][mg.core_nodes],
... np.array([ 0.60653066, 0.36787944, 0.22313016]))
True
>>> DDdiff.run_one_step(2.)
>>> np.allclose(
... mg.at_node['topographic__elevation'][mg.core_nodes],
... np.array([ 1.47730244, 2.28949856, 3.17558975]))
True
>>> np.allclose(mg.at_node['bedrock__elevation'][mg.core_nodes],
... np.array([-0.71306132, 0.26424112, 1.05373968]))
True
>>> np.allclose(mg.at_node['soil__depth'], z - BRz)
True
Now, we’ll test that changing the transport decay depth behaves as expected.
>>> mg = RasterModelGrid((3, 5))
>>> soilTh = mg.add_zeros("soil__depth", at="node")
>>> z = mg.add_zeros("topographic__elevation", at="node")
>>> BRz = mg.add_zeros("bedrock__elevation", at="node")
>>> z += mg.node_x.copy()**0.5
>>> BRz = z.copy() - 1.0
>>> soilTh[:] = z - BRz
>>> expweath = ExponentialWeatherer(mg)
>>> DDdiff = DepthDependentDiffuser(mg, soil_transport_decay_depth = 0.1)
>>> DDdiff.run_one_step(1)
>>> soil_decay_depth_point1 = mg.at_node['topographic__elevation'][mg.core_nodes]
>>> z[:] = 0
>>> z += mg.node_x.copy()**0.5
>>> BRz = z.copy() - 1.0
>>> soilTh[:] = z - BRz
>>> DDdiff = DepthDependentDiffuser(mg, soil_transport_decay_depth = 1.0)
>>> DDdiff.run_one_step(1)
>>> soil_decay_depth_1 = mg.at_node['topographic__elevation'][mg.core_nodes]
>>> np.greater(soil_decay_depth_1[1], soil_decay_depth_point1[1])
False
References
Required Software Citation(s) Specific to this Component
Barnhart, K., Glade, R., Shobe, C., Tucker, G. (2019). Terrainbento 1.0: a Python package for multi-model analysis in long-term drainage basin evolution. Geoscientific Model Development 12(4), 1267–1297. https://dx.doi.org/10.5194/gmd-12-1267-2019
Additional References
Johnstone, S., Hilley, G. (2015). Lithologic control on the form of soil-mantled hillslopes Geology 43(1), 83-86. https://dx.doi.org/10.1130/g36052.1
| Parameters: |
|
|---|
DepthDependentTaylorDiffuser(grid, linear_diffusivity=1.0, slope_crit=1.0, soil_transport_decay_depth=1.0, nterms=2, dynamic_dt=False, if_unstable='pass', courant_factor=0.2)[source]¶Bases: landlab.core.model_component.Component
This component implements a depth-dependent Taylor series diffusion rule, combining concepts of Ganti et al. (2012) and Johnstone and Hilley (2014).
Hillslope sediment flux uses a Taylor Series expansion of the Andrews- Bucknam formulation of nonlinear hillslope flux derived following following Ganti et al., 2012 with a depth dependent component inspired Johnstone and Hilley (2014). The flux \(q_s\) is given as:
where \(D\) is is the diffusivity, \(S\) is the slope (defined as negative downward), \(S_c\) is the critical slope, \(n\) is the number of terms, \(H\) is the soil depth on links, and \(H^*\) is the soil transport decay depth.
The default behavior uses two terms to produce a slope dependence as described by Equation 6 of Ganti et al., (2012).
This component will ignore soil thickness located at non-core nodes.
Examples
First lets make a simple example with flat topography.
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import ExponentialWeatherer
>>> from landlab.components import DepthDependentTaylorDiffuser
>>> mg = RasterModelGrid((5, 5))
>>> soilTh = mg.add_zeros('node', 'soil__depth')
>>> z = mg.add_zeros('node', 'topographic__elevation')
>>> BRz = mg.add_zeros('node', 'bedrock__elevation')
>>> expweath = ExponentialWeatherer(mg)
>>> DDdiff = DepthDependentTaylorDiffuser(mg)
>>> expweath.calc_soil_prod_rate()
>>> np.allclose(mg.at_node['soil_production__rate'][mg.core_nodes], 1.)
True
>>> DDdiff.run_one_step(2.)
>>> np.allclose(mg.at_node['topographic__elevation'][mg.core_nodes], 0.)
True
>>> np.allclose(mg.at_node['bedrock__elevation'][mg.core_nodes], -2.)
True
>>> np.allclose(mg.at_node['soil__depth'][mg.core_nodes], 2.)
True
Now a more complicated example with a slope.
>>> mg = RasterModelGrid((3, 5))
>>> soilTh = mg.add_zeros('node', 'soil__depth')
>>> z = mg.add_zeros('node', 'topographic__elevation')
>>> BRz = mg.add_zeros('node', 'bedrock__elevation')
>>> z += mg.node_x.copy()
>>> BRz += mg.node_x / 2.
>>> soilTh[:] = z - BRz
>>> expweath = ExponentialWeatherer(mg)
>>> DDdiff = DepthDependentTaylorDiffuser(mg)
>>> expweath.calc_soil_prod_rate()
>>> np.allclose(
... mg.at_node['soil_production__rate'][mg.core_nodes],
... np.array([ 0.60653066, 0.36787944, 0.22313016]))
True
>>> DDdiff.run_one_step(0.1)
>>> np.allclose(
... mg.at_node['topographic__elevation'][mg.core_nodes],
... np.array([ 1.04773024, 2.02894986, 3.01755898]))
True
>>> np.allclose(mg.at_node['bedrock__elevation'][mg.core_nodes],
... np.array([ 0.43934693, 0.96321206, 1.47768698]))
True
>>> np.allclose(mg.at_node['soil__depth'], z - BRz)
True
The DepthDependentTaylorDiffuser makes and moves soil at a rate proportional to slope, this means that there is a characteristic time scale for soil transport and an associated stability criteria for the timestep. The maximum characteristic time scale, \(De_{max}\), is given as a function of the hillslope diffustivity, \(D\), the maximum slope, \(S_{max}\), and the critical slope \(S_c\).
The maximum stable time step is given by
Where the courant factor is a user defined scale (default is 0.2)
The DepthDependentTaylorDiffuser has a boolean flag that permits a user to be warned if timesteps are too large for the slopes in the model grid (if_unstable = ‘warn’) and a boolean flag that turns on dynamic timestepping (dynamic_dt = False).
>>> DDdiff = DepthDependentTaylorDiffuser(mg, if_unstable='warn')
>>> DDdiff.run_one_step(2.)
Topographic slopes are high enough such that the Courant condition is
exceeded AND you have not selected dynamic timestepping with
dynamic_dt=True. This may lead to infinite and/or nan values for slope,
elevation, and soil depth. Consider using a smaller time step or dynamic
timestepping. The Courant condition recommends a timestep of
0.0953407607307 or smaller.
Alternatively you can specify if_unstable=’raise’, and a Runtime Error will be raised if this condition is not met.
Next, lets do an example with dynamic timestepping.
>>> mg = RasterModelGrid((3, 5))
>>> soilTh = mg.add_zeros('node', 'soil__depth')
>>> z = mg.add_zeros('node', 'topographic__elevation')
>>> BRz = mg.add_zeros('node', 'bedrock__elevation')
We’ll use a steep slope and very little soil.
>>> z += mg.node_x.copy()**2
>>> BRz = z.copy() - 1.0
>>> soilTh[:] = z - BRz
>>> expweath = ExponentialWeatherer(mg)
Lets try to move the soil with a large timestep. Without dynamic time steps, this gives a warning that we’ve exceeded the dynamic timestep size and should use a smaller timestep. We could either use the smaller timestep, or specify that we want to use a dynamic timestep.
>>> DDdiff = DepthDependentTaylorDiffuser(mg, if_unstable='warn', dynamic_dt=False)
>>> expweath.calc_soil_prod_rate()
>>> DDdiff.run_one_step(10)
Topographic slopes are high enough such that the Courant condition is
exceeded AND you have not selected dynamic timestepping with
dynamic_dt=True. This may lead to infinite and/or nan values for slope,
elevation, and soil depth. Consider using a smaller time step or dynamic
timestepping. The Courant condition recommends a timestep of
0.004 or smaller.
Now, we’ll re-build the grid and do the same example with dynamic timesteps.
>>> mg = RasterModelGrid((3, 5))
>>> soilTh = mg.add_zeros('node', 'soil__depth')
>>> z = mg.add_zeros('node', 'topographic__elevation')
>>> BRz = mg.add_zeros('node', 'bedrock__elevation')
>>> z += mg.node_x.copy()**2
>>> BRz = z.copy() - 1.0
>>> soilTh[:] = z - BRz
>>> expweath = ExponentialWeatherer(mg)
>>> DDdiff = DepthDependentTaylorDiffuser(mg, if_unstable='warn', dynamic_dt=True)
>>> expweath.calc_soil_prod_rate()
>>> DDdiff.run_one_step(10)
>>> np.any(np.isnan(z))
False
Now, we’ll test that changing the transport decay depth behaves as expected.
>>> mg = RasterModelGrid((3, 5))
>>> soilTh = mg.add_zeros('node', 'soil__depth')
>>> z = mg.add_zeros('node', 'topographic__elevation')
>>> BRz = mg.add_zeros('node', 'bedrock__elevation')
>>> z += mg.node_x.copy()**0.5
>>> BRz = z.copy() - 1.0
>>> soilTh[:] = z - BRz
>>> expweath = ExponentialWeatherer(mg)
>>> DDdiff = DepthDependentTaylorDiffuser(mg, soil_transport_decay_depth = 0.1)
>>> DDdiff.run_one_step(1)
>>> soil_decay_depth_point1 = mg.at_node['topographic__elevation'][mg.core_nodes]
>>> z[:] = 0
>>> z += mg.node_x.copy()**0.5
>>> BRz = z.copy() - 1.0
>>> soilTh[:] = z - BRz
>>> DDdiff = DepthDependentTaylorDiffuser(mg, soil_transport_decay_depth = 1.0)
>>> DDdiff.run_one_step(1)
>>> soil_decay_depth_1 = mg.at_node['topographic__elevation'][mg.core_nodes]
>>> np.greater(soil_decay_depth_1[1], soil_decay_depth_point1[1])
False
References
Required Software Citation(s) Specific to this Component
Barnhart, K., Glade, R., Shobe, C., Tucker, G. (2019). Terrainbento 1.0: a Python package for multi-model analysis in long-term drainage basin evolution. Geoscientific Model Development 12(4), 1267–1297. https://dx.doi.org/10.5194/gmd-12-1267-2019
Additional References
Ganti, V., Passalacqua, P., Foufoula-Georgiou, E. (2012). A sub-grid scale closure for nonlinear hillslope sediment transport models Journal of Geophysical Research: Earth Surface 117(F2). https://dx.doi.org/10.1029/2011jf002181
Johnstone, S., Hilley, G. (2015). Lithologic control on the form of soil-mantled hillslopes Geology 43(1), 83-86. https://dx.doi.org/10.1130/g36052.1
Initialize the DepthDependentTaylorDiffuser.
| Parameters: |
|
|---|
__init__(grid, linear_diffusivity=1.0, slope_crit=1.0, soil_transport_decay_depth=1.0, nterms=2, dynamic_dt=False, if_unstable='pass', courant_factor=0.2)[source]¶Initialize the DepthDependentTaylorDiffuser.
| Parameters: |
|
|---|
DepthSlopeProductErosion(grid, k_e=0.001, fluid_density=1000.0, g=9.80665, a_exp=1.0, tau_crit=0.0, uplift_rate=0.0, slope='topographic__slope')[source]¶Bases: landlab.core.model_component.Component
Landlab component that simulates detachment-limited river erosion.
This component calculates changes in elevation in response to vertical incision.
Calculate detachment limited erosion rate on nodes using the shear stress equation, solved using the depth slope product.
Landlab component that generalizes the detachment limited erosion equation, primarily to be coupled to the the Landlab OverlandFlow component.
This component adjusts topographic elevation and is contained in the landlab.components.detachment_ltd_erosion folder.
| Parameters: |
|
|---|
__init__(grid, k_e=0.001, fluid_density=1000.0, g=9.80665, a_exp=1.0, tau_crit=0.0, uplift_rate=0.0, slope='topographic__slope')[source]¶Calculate detachment limited erosion rate on nodes using the shear stress equation, solved using the depth slope product.
Landlab component that generalizes the detachment limited erosion equation, primarily to be coupled to the the Landlab OverlandFlow component.
This component adjusts topographic elevation and is contained in the landlab.components.detachment_ltd_erosion folder.
| Parameters: |
|
|---|
dz¶Magnitude of change of the topographic__elevation due to erosion [L].
DetachmentLtdErosion(grid, K_sp=2e-05, m_sp=0.5, n_sp=1.0, uplift_rate=0.0, entrainment_threshold=0.0, slope='topographic__slope')[source]¶Bases: landlab.core.model_component.Component
Landlab component that simulates detachment-limited river erosion.
This component calculates changes in elevation in response to vertical incision.
Calculate detachment limited erosion rate on nodes.
Landlab component that generalizes the detachment limited erosion equation, primarily to be coupled to the the Landlab OverlandFlow component.
This component adjusts topographic elevation.
| Parameters: |
|
|---|
__init__(grid, K_sp=2e-05, m_sp=0.5, n_sp=1.0, uplift_rate=0.0, entrainment_threshold=0.0, slope='topographic__slope')[source]¶Calculate detachment limited erosion rate on nodes.
Landlab component that generalizes the detachment limited erosion equation, primarily to be coupled to the the Landlab OverlandFlow component.
This component adjusts topographic elevation.
| Parameters: |
|
|---|
DischargeDiffuser(grid, slope=0.25, flat_thresh=0.0001)[source]¶Bases: landlab.core.model_component.Component
Diffuse sediment proportional to an implicit water discharge value.
This class implements Voller, Hobley, and Paola’s scheme for sediment diffusion, where the diffusivity of the sediment is proportional to the local discharge of water. The method works by solving for a potential field describing the water discharge at all nodes on the grid, which enforces both mass conservation and flow downhill along topographic gradients. This routine is designed to construct sediment fans.
Note that both the water and sediment discharges are calculated together within the component.
The algorithm uses a rule that looks like:
q_sed = q_water * (S - S_crit)
where S_crit is a critical slope threshold. [MODIFY THIS]
It is VITAL you initialize this component AFTER setting boundary conditions.
The primary method of this class is run_one_step.
Notes
This is a “research grade” component, and is subject to dramatic change with little warning. No guarantees are made regarding its accuracy or utility. It is not recommended for user use yet!
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
None Listed
| Parameters: | grid (ModelGrid) – A grid. |
|---|
discharges_at_links¶Return the discharges at links.
Note that if diagonal routing, this will return number_of_d8_links. Otherwise, it will be number_of_links.
DrainageDensity(grid, channel__mask=None, area_coefficient=None, slope_coefficient=None, area_exponent=None, slope_exponent=None, channelization_threshold=None)[source]¶Bases: landlab.core.model_component.Component
Calculate drainage density over a DEM.
Landlab component that implements the distance to channel algorithm of Tucker et al., 2001.
This component requires EITHER a channel__mask array with 1’s
where channels exist and 0’s elsewhere, OR a set of coefficients
and exponents for a slope-area relationship and a
channelization threshold to compare against that relationship.
If an array is provided it MUST be of type np.uint8. See the example
below for how to make such an array.
The channel__mask array will be assigned to an at-node field with the
name channel__mask. If the channel__mask was originaly created from a
passed array, a user can update this array to change the mask.
If the channel__mask is created using an area coefficent,
slope coefficient, area exponent, slope exponent, and channelization
threshold, the location of the mask will be re-update when
calculate_drainage_density is called.
If an area coefficient, \(C_A\), a slope coefficent, \(C_S\), an area exponent, \(m_r\), a slope exponent, \(n_r\), and channelization threshold \(T_C\) are provided, nodes that meet the criteria
where \(A\) is the drainage density and \(S\) is the local slope, will be marked as channel nodes.
The calculate_drainage_density function returns drainage density for the
model domain. This function calculates the distance from every node to the
nearest channel node \(L\) along the flow line of steepest descent
(assuming D8 routing).
This component stores this distance a field, called:
surface_to_channel__minimum_distance. The drainage density is then
calculated (after Tucker et al., 2001):
where \(\overline{L}\) is the mean L for the model domain.
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator, FastscapeEroder
>>> mg = RasterModelGrid((10, 10))
>>> _ = mg.add_zeros('node', 'topographic__elevation')
>>> np.random.seed(50)
>>> noise = np.random.rand(100)
>>> mg.at_node['topographic__elevation'] += noise
>>> mg.at_node['topographic__elevation'] # doctest: +NORMALIZE_WHITESPACE
array([ 0.49460165, 0.2280831 , 0.25547392, 0.39632991, 0.3773151 ,
0.99657423, 0.4081972 , 0.77189399, 0.76053669, 0.31000935,
0.3465412 , 0.35176482, 0.14546686, 0.97266468, 0.90917844,
0.5599571 , 0.31359075, 0.88820004, 0.67457307, 0.39108745,
0.50718412, 0.5241035 , 0.92800093, 0.57137307, 0.66833757,
0.05225869, 0.3270573 , 0.05640164, 0.17982769, 0.92593317,
0.93801522, 0.71409271, 0.73268761, 0.46174768, 0.93132927,
0.40642024, 0.68320577, 0.64991587, 0.59876518, 0.22203939,
0.68235717, 0.8780563 , 0.79671726, 0.43200225, 0.91787822,
0.78183368, 0.72575028, 0.12485469, 0.91630845, 0.38771099,
0.29492955, 0.61673141, 0.46784623, 0.25533891, 0.83899589,
0.1786192 , 0.22711417, 0.65987645, 0.47911625, 0.07344734,
0.13896007, 0.11230718, 0.47778497, 0.54029623, 0.95807105,
0.58379231, 0.52666409, 0.92226269, 0.91925702, 0.25200886,
0.68263261, 0.96427612, 0.22696165, 0.7160172 , 0.79776011,
0.9367512 , 0.8537225 , 0.42154581, 0.00543987, 0.03486533,
0.01390537, 0.58890993, 0.3829931 , 0.11481895, 0.86445401,
0.82165703, 0.73749168, 0.84034417, 0.4015291 , 0.74862 ,
0.55962945, 0.61323757, 0.29810165, 0.60237917, 0.42567684,
0.53854438, 0.48672986, 0.49989164, 0.91745948, 0.26287702])
>>> fr = FlowAccumulator(mg, flow_director='D8')
>>> fsc = FastscapeEroder(mg, K_sp=.01, m_sp=.5, n_sp=1)
>>> for x in range(100):
... fr.run_one_step()
... fsc.run_one_step(dt = 10.0)
... mg.at_node['topographic__elevation'][mg.core_nodes] += .01
>>> channels = np.array(mg.at_node['drainage_area'] > 5, dtype=np.uint8)
>>> dd = DrainageDensity(mg, channel__mask=channels)
>>> mean_drainage_density = dd.calculate_drainage_density()
>>> np.isclose(mean_drainage_density, 0.3831100571)
True
Alternatively you can pass a set of coefficients to identify the channel mask. Next shows the same example as above, but with these coefficients provided.
>>> mg = RasterModelGrid((10, 10))
>>> _ = mg.add_zeros('node', 'topographic__elevation')
>>> np.random.seed(50)
>>> noise = np.random.rand(100)
>>> mg.at_node['topographic__elevation'] += noise
>>> fr = FlowAccumulator(mg, flow_director='D8')
>>> fsc = FastscapeEroder(mg, K_sp=.01, m_sp=.5, n_sp=1)
>>> for x in range(100):
... fr.run_one_step()
... fsc.run_one_step(dt = 10.0)
... mg.at_node['topographic__elevation'][mg.core_nodes] += .01
>>> channels = np.array(mg.at_node['drainage_area'] > 5, dtype=np.uint8)
>>> dd = DrainageDensity(mg,
... area_coefficient=1.0,
... slope_coefficient=1.0,
... area_exponent=1.0,
... slope_exponent=0.0,
... channelization_threshold=5)
>>> mean_drainage_density = dd.calculate_drainage_density()
>>> np.isclose(mean_drainage_density, 0.3831100571)
True
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
Tucker, G., Catani, F., Rinaldo, A., Bras, R. (2001). Statistical analysis of drainage density from digital terrain data. Geomorphology 36(3-4), 187-202. https://dx.doi.org/10.1016/s0169-555x(00)00056-8
Initialize the DrainageDensity component.
| Parameters: |
|
|---|
__init__(grid, channel__mask=None, area_coefficient=None, slope_coefficient=None, area_exponent=None, slope_exponent=None, channelization_threshold=None)[source]¶Initialize the DrainageDensity component.
| Parameters: |
|
|---|
calculate_drainage_density()[source]¶Calculate drainage density.
If the channel mask is defined based on slope and area coefficients, it will be update based on the current drainage area and slope fields.
| Returns: | landscape_drainage_density – Drainage density over the model domain. |
|---|---|
| Return type: | float (1/m) |
ErosionDeposition(grid, K=0.002, v_s=1.0, m_sp=0.5, n_sp=1.0, sp_crit=0.0, F_f=0.0, discharge_field='surface_water__discharge', solver='basic', dt_min=0.001, **kwds)[source]¶Bases: landlab.components.erosion_deposition.generalized_erosion_deposition._GeneralizedErosionDeposition
Erosion-Deposition model in the style of Davy and Lague (2009). It uses a mass balance approach across the total sediment mass both in the bed and in transport coupled with explicit representation of the sediment transport lengthscale (the “xi-q” model) to derive a range of erosional and depositional responses in river channels.
This implementation is close to the Davy & Lague scheme, with a few deviations:
- A fraction of the eroded sediment is permitted to enter the wash load, and lost to the mass balance (F_f).
- Here an incision threshold \(\omega\) is permitted, where it was not by Davy & Lague. It is implemented with an exponentially smoothed form to prevent discontinuities in the parameter space. See the
StreamPowerSmoothThresholdEroderfor more documentation.- This component uses an “effective” settling velocity, v_s, as one of its inputs. This parameter is simply equal to Davy & Lague’s d_star * V dimensionless number.
Erosion of the bed follows a stream power formulation, i.e.,
Note that the transition between transport-limited and detachment-limited behavior is controlled by the dimensionless ratio (v_s/r) where r is the runoff ratio (Q=Ar). r can be changed in the flow accumulation component but is not changed within ErosionDeposition. Because the runoff ratio r is not changed within the ErosionDeposition component, v_s becomes the parameter that fundamentally controls response style. Very small v_s will lead to a detachment-limited response style, very large v_s will lead to a transport-limited response style. v_s == 1 means equal contributions from transport and erosion, and a hybrid response as described by Davy & Lague.
Unlike other some other fluvial erosion componets in Landlab, in this
component (and SPACE) no erosion occurs
in depressions or in areas with adverse slopes. There is no ability to
pass a keyword argument erode_flooded_nodes.
If a depressions are handled (as indicated by the presence of the field “flood_status_code” at nodes), then deposition occurs throughout the depression and sediment is passed out of the depression. Where pits are encountered, then all sediment is deposited at that node only.
A note about sediment porosity: Prior to Landlab v2.0 this component took a
porositiy keyworkd argument phi. For an explaination of why it no
longer does (including a mathematical derivation), see
Pull Request 1186.
If phi is passed to this component a value error will be raised.
Component written by C. Shobe, K. Barnhart, and G. Tucker.
References
Required Software Citation(s) Specific to this Component
Barnhart, K., Glade, R., Shobe, C., Tucker, G. (2019). Terrainbento 1.0: a Python package for multi-model analysis in long-term drainage basin evolution. Geoscientific Model Development 12(4), 1267–1297. https://dx.doi.org/10.5194/gmd-12-1267-2019
Additional References
Davy, P., Lague, D. (2009). Fluvial erosion/transport equation of landscape evolution models revisited Journal of Geophysical Research 114(F3), F03007. https://dx.doi.org/10.1029/2008jf001146
Initialize the ErosionDeposition model.
| Parameters: |
|
|---|
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> from landlab.components import DepressionFinderAndRouter
>>> from landlab.components import ErosionDeposition
>>> from landlab.components import FastscapeEroder
>>> np.random.seed(seed = 5000)
>>> nr = 5
>>> nc = 5
>>> dx = 10
>>> mg = RasterModelGrid((nr, nc), xy_spacing=10.0)
>>> _ = mg.add_zeros('node', 'topographic__elevation')
>>> mg['node']['topographic__elevation'] += (mg.node_y/10 +
... mg.node_x/10 + np.random.rand(len(mg.node_y)) / 10)
>>> mg.set_closed_boundaries_at_grid_edges(bottom_is_closed=True,
... left_is_closed=True,
... right_is_closed=True,
... top_is_closed=True)
>>> mg.set_watershed_boundary_condition_outlet_id(0, mg['node']['topographic__elevation'], -9999.)
>>> fsc_dt = 100.
>>> ed_dt = 1.
Check initial topography
>>> mg.at_node['topographic__elevation'] # doctest: +NORMALIZE_WHITESPACE
array([ 0.02290479, 1.03606698, 2.0727653 , 3.01126678, 4.06077707,
1.08157495, 2.09812694, 3.00637448, 4.07999597, 5.00969486,
2.04008677, 3.06621577, 4.09655859, 5.04809001, 6.02641123,
3.05874171, 4.00585786, 5.0595697 , 6.04425233, 7.05334077,
4.05922478, 5.0409473 , 6.07035008, 7.0038935 , 8.01034357])
Instantiate Fastscape eroder, flow router, and depression finder
>>> fr = FlowAccumulator(mg, flow_director='D8')
>>> df = DepressionFinderAndRouter(mg)
>>> fsc = FastscapeEroder(
... mg,
... K_sp=.001,
... m_sp=.5,
... n_sp=1)
Burn in an initial drainage network using the Fastscape eroder:
>>> for x in range(100):
... fr.run_one_step()
... df.map_depressions()
... flooded = np.where(df.flood_status==3)[0]
... fsc.run_one_step(dt = fsc_dt)
... mg.at_node['topographic__elevation'][0] -= 0.001 #uplift
Instantiate the E/D component:
>>> ed = ErosionDeposition(
... mg,
... K=0.00001,
... v_s=0.001,
... m_sp=0.5,
... n_sp = 1.0,
... sp_crit=0)
Now run the E/D component for 2000 short timesteps:
>>> for x in range(2000): #E/D component loop
... fr.run_one_step()
... df.map_depressions()
... ed.run_one_step(dt = ed_dt)
... mg.at_node['topographic__elevation'][0] -= 2e-4 * ed_dt
Now we test to see if topography is right:
>>> np.around(mg.at_node['topographic__elevation'], decimals=3) # doctest: +NORMALIZE_WHITESPACE
array([-0.477, 1.036, 2.073, 3.011, 4.061, 1.082, -0.08 , -0.065,
-0.054, 5.01 , 2.04 , -0.065, -0.065, -0.053, 6.026, 3.059,
-0.054, -0.053, -0.035, 7.053, 4.059, 5.041, 6.07 , 7.004,
8.01 ])
K¶Erodibility (units depend on m_sp).
__init__(grid, K=0.002, v_s=1.0, m_sp=0.5, n_sp=1.0, sp_crit=0.0, F_f=0.0, discharge_field='surface_water__discharge', solver='basic', dt_min=0.001, **kwds)[source]¶Initialize the ErosionDeposition model.
| Parameters: |
|
|---|
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> from landlab.components import DepressionFinderAndRouter
>>> from landlab.components import ErosionDeposition
>>> from landlab.components import FastscapeEroder
>>> np.random.seed(seed = 5000)
>>> nr = 5
>>> nc = 5
>>> dx = 10
>>> mg = RasterModelGrid((nr, nc), xy_spacing=10.0)
>>> _ = mg.add_zeros('node', 'topographic__elevation')
>>> mg['node']['topographic__elevation'] += (mg.node_y/10 +
... mg.node_x/10 + np.random.rand(len(mg.node_y)) / 10)
>>> mg.set_closed_boundaries_at_grid_edges(bottom_is_closed=True,
... left_is_closed=True,
... right_is_closed=True,
... top_is_closed=True)
>>> mg.set_watershed_boundary_condition_outlet_id(0, mg['node']['topographic__elevation'], -9999.)
>>> fsc_dt = 100.
>>> ed_dt = 1.
Check initial topography
>>> mg.at_node['topographic__elevation'] # doctest: +NORMALIZE_WHITESPACE
array([ 0.02290479, 1.03606698, 2.0727653 , 3.01126678, 4.06077707,
1.08157495, 2.09812694, 3.00637448, 4.07999597, 5.00969486,
2.04008677, 3.06621577, 4.09655859, 5.04809001, 6.02641123,
3.05874171, 4.00585786, 5.0595697 , 6.04425233, 7.05334077,
4.05922478, 5.0409473 , 6.07035008, 7.0038935 , 8.01034357])
Instantiate Fastscape eroder, flow router, and depression finder
>>> fr = FlowAccumulator(mg, flow_director='D8')
>>> df = DepressionFinderAndRouter(mg)
>>> fsc = FastscapeEroder(
... mg,
... K_sp=.001,
... m_sp=.5,
... n_sp=1)
Burn in an initial drainage network using the Fastscape eroder:
>>> for x in range(100):
... fr.run_one_step()
... df.map_depressions()
... flooded = np.where(df.flood_status==3)[0]
... fsc.run_one_step(dt = fsc_dt)
... mg.at_node['topographic__elevation'][0] -= 0.001 #uplift
Instantiate the E/D component:
>>> ed = ErosionDeposition(
... mg,
... K=0.00001,
... v_s=0.001,
... m_sp=0.5,
... n_sp = 1.0,
... sp_crit=0)
Now run the E/D component for 2000 short timesteps:
>>> for x in range(2000): #E/D component loop
... fr.run_one_step()
... df.map_depressions()
... ed.run_one_step(dt = ed_dt)
... mg.at_node['topographic__elevation'][0] -= 2e-4 * ed_dt
Now we test to see if topography is right:
>>> np.around(mg.at_node['topographic__elevation'], decimals=3) # doctest: +NORMALIZE_WHITESPACE
array([-0.477, 1.036, 2.073, 3.011, 4.061, 1.082, -0.08 , -0.065,
-0.054, 5.01 , 2.04 , -0.065, -0.065, -0.053, 6.026, 3.059,
-0.054, -0.053, -0.035, 7.053, 4.059, 5.041, 6.07 , 7.004,
8.01 ])
ExponentialWeatherer(grid, soil_production__maximum_rate=1.0, soil_production__decay_depth=1.0)[source]¶Bases: landlab.core.model_component.Component
This component implements exponential weathering of bedrock on hillslopes. Uses exponential soil production function in the style of Ahnert (1976).
Consider that \(w_0\) is the maximum soil production rate and that \(d^*\) is the characteristic soil production depth. The soil production rate \(w\) is given as a function of the soil depth \(d\),
The ExponentialWeatherer only calculates soil production at core nodes.
An alternative version which uses the analytical integral of
production through time is available at the component
ExponentialWeathererIntegrated.
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import ExponentialWeatherer
>>> mg = RasterModelGrid((5, 5))
>>> soilz = mg.add_zeros("soil__depth", at="node")
>>> soilrate = mg.add_ones("soil_production__rate", at="node")
>>> expw = ExponentialWeatherer(mg)
>>> expw.calc_soil_prod_rate()
>>> np.allclose(mg.at_node['soil_production__rate'], 1.)
True
References
Required Software Citation(s) Specific to this Component
Barnhart, K., Glade, R., Shobe, C., Tucker, G. (2019). Terrainbento 1.0: a Python package for multi-model analysis in long-term drainage basin evolution. Geoscientific Model Development 12(4), 1267–1297. https://dx.doi.org/10.5194/gmd-12-1267-2019
Additional References
Ahnert, F. (1976). Brief description of a comprehensive three-dimensional process-response model of landform development Z. Geomorphol. Suppl. 25, 29 - 49.
Armstrong, A. (1976). A three dimensional simulation of slope forms. Zeitschrift für Geomorphologie 25, 20 - 28.
| Parameters: |
|
|---|
__init__(grid, soil_production__maximum_rate=1.0, soil_production__decay_depth=1.0)[source]¶| Parameters: |
|
|---|
maximum_weathering_rate¶Maximum rate of weathering (m/yr).
ExponentialWeathererIntegrated(grid, soil_production__maximum_rate=1.0, soil_production__decay_depth=1.0, soil_production__expansion_factor=1.0)[source]¶Bases: landlab.core.model_component.Component
This component implements exponential weathering of bedrock on hillslopes. Uses exponential soil production function in the style of Ahnert (1976).
Consider that \(w_0\) is the maximum soil production rate and that \(d^*\) is the characteristic soil production depth. The soil production rate \(w\) is given as a function of the soil depth \(d\),
The ExponentialWeathererIntegrated uses the analytical solution for the amount of soil produced by an exponential weathering function over a timestep dt, and returns both the thickness of bedrock weathered and the thickness of soil produced. The solution accounts for the reduction in rate over the timestep due to the increasing depth. This enables accuracy over arbitrarily large timesteps, and better compatiblity with the run_one_step() interface.
Compared to ‘ExponentialWeatherer’, upon which it is based…
This is left as an exercise for the model driver, as different applications may want to integrate soil depth and weathering in different sequences among other processes.
ExponentialWeatherer,
just import and instantiate this one instead and existing code
should work with no side effects other than the creation of the
two additional (zeros) output fields.Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import ExponentialWeathererIntegrated
>>> mg = RasterModelGrid((5, 5))
>>> soilz = mg.add_zeros("soil__depth", at="node")
>>> soilrate = mg.add_ones("soil_production__rate", at="node")
>>> expw = ExponentialWeathererIntegrated(mg)
>>> dt = 1000
>>> expw.run_one_step(dt)
>>> np.allclose(mg.at_node['soil_production__rate'][mg.core_nodes], 1.)
True
>>> np.allclose(mg.at_node['soil_production__dt_produced_depth'][mg.core_nodes], 6.9088)
True
References
Required Software Citation(s) Specific to this Component
Barnhart, K., Glade, R., Shobe, C., Tucker, G. (2019). Terrainbento 1.0: a Python package for multi-model analysis in long-term drainage basin evolution. Geoscientific Model Development 12(4), 1267–1297. https://dx.doi.org/10.5194/gmd-12-1267-2019
Additional References
Ahnert, F. (1976). Brief description of a comprehensive three-dimensional process-response model of landform development Z. Geomorphol. Suppl. 25, 29 - 49.
Armstrong, A. (1976). A three dimensional simulation of slope forms. Zeitschrift für Geomorphologie 25, 20 - 28.
| Parameters: |
|
|---|
__init__(grid, soil_production__maximum_rate=1.0, soil_production__decay_depth=1.0, soil_production__expansion_factor=1.0)[source]¶| Parameters: |
|
|---|
maximum_weathering_rate¶Maximum rate of weathering (m/yr).
FastscapeEroder(grid, K_sp=0.001, m_sp=0.5, n_sp=1.0, threshold_sp=0.0, discharge_field='drainage_area', erode_flooded_nodes=True)[source]¶Bases: landlab.core.model_component.Component
Fastscape stream power erosion.
This class uses the Braun-Willett Fastscape approach to calculate the amount of erosion at each node in a grid, following a stream power framework. This should allow it to be stable against larger timesteps than an explicit stream power scheme.
Note that although this scheme is nominally implicit, and will reach a numerically-correct solution under topographic steady state regardless of timestep length, the accuracy of transient solutions is not timestep independent (see Braun & Willett 2013, Appendix B for further details). Although the scheme remains significantly more robust and permits longer timesteps than a traditional explicit solver under such conditions, it is still possible to create numerical instability through use of too long a timestep while using this component. The user is cautioned to check their implementation is behaving stably before fully trusting it.
Stream power erosion is implemented as:
if \(K A ^ m S ^ n > \textit{threshold_sp}\), and:
if \(K A^m S^n <= \textit{threshold_sp}\).
This module assumes you have already run
landlab.components.flow_accum.flow_accumulator.FlowAccumulator.run_one_step
in the same timestep. It looks for ‘flow__upstream_node_order’,
‘flow__link_to_receiver_node’, ‘drainage_area’, ‘flow__receiver_node’, and
‘topographic__elevation’ at the nodes in the grid. ‘drainage_area’ should
be in area upstream, not volume (i.e., set runoff_rate=1.0 when calling
FlowAccumulator.run_one_step).
The primary method of this class is run_one_step.
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator, FastscapeEroder
>>> grid = RasterModelGrid((5, 5), xy_spacing=10.)
>>> z = np.array([7., 7., 7., 7., 7.,
... 7., 5., 3.2, 6., 7.,
... 7., 2., 3., 5., 7.,
... 7., 1., 1.9, 4., 7.,
... 7., 0., 7., 7., 7.])
>>> z = grid.add_field('topographic__elevation', z, at='node')
>>> fr = FlowAccumulator(grid, flow_director='D8')
>>> sp = FastscapeEroder(grid, K_sp=1.)
>>> fr.run_one_step()
>>> sp.run_one_step(dt=1.)
>>> z # doctest: +NORMALIZE_WHITESPACE
array([ 7. , 7. , 7. , 7. , 7. ,
7. , 2.92996598, 2.02996598, 4.01498299, 7. ,
7. , 0.85993197, 1.87743897, 3.28268321, 7. ,
7. , 0.28989795, 0.85403051, 2.42701526, 7. ,
7. , 0. , 7. , 7. , 7. ])
>>> grid = RasterModelGrid((3, 7), xy_spacing=1.)
>>> z = np.array(grid.node_x ** 2.)
>>> z = grid.add_field('topographic__elevation', z, at='node')
>>> grid.status_at_node[grid.nodes_at_left_edge] = grid.BC_NODE_IS_FIXED_VALUE
>>> grid.status_at_node[grid.nodes_at_top_edge] = grid.BC_NODE_IS_CLOSED
>>> grid.status_at_node[grid.nodes_at_bottom_edge] = grid.BC_NODE_IS_CLOSED
>>> grid.status_at_node[grid.nodes_at_right_edge] = grid.BC_NODE_IS_CLOSED
>>> fr = FlowAccumulator(grid, flow_director='D8')
>>> sp = FastscapeEroder(grid, K_sp=0.1, m_sp=0., n_sp=2.,
... threshold_sp=2.)
>>> fr.run_one_step()
>>> sp.run_one_step(dt=10.)
>>> z.reshape(grid.shape)[1, :] # doctest: +NORMALIZE_WHITESPACE
array([ 0. , 1. , 4. , 8.52493781,
13.29039716, 18.44367965, 36. ])
>>> grid = RasterModelGrid((3, 7), xy_spacing=1.)
>>> z = np.array(grid.node_x ** 2.)
>>> z = grid.add_field('topographic__elevation', z, at='node')
>>> grid.status_at_node[grid.nodes_at_left_edge] = grid.BC_NODE_IS_FIXED_VALUE
>>> grid.status_at_node[grid.nodes_at_top_edge] = grid.BC_NODE_IS_CLOSED
>>> grid.status_at_node[grid.nodes_at_bottom_edge] = grid.BC_NODE_IS_CLOSED
>>> grid.status_at_node[grid.nodes_at_right_edge] = grid.BC_NODE_IS_CLOSED
>>> cell_area = 1.0
>>> fr = FlowAccumulator(grid, flow_director='D8', runoff_rate=2.0)
>>> grid.at_node["water__unit_flux_in"]
array([ 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2.,
2., 2., 2., 2., 2., 2., 2., 2.])
>>> K_field = grid.ones(at='node') # K can be a field
>>> sp = FastscapeEroder(grid,
... K_sp=K_field,
... m_sp=1.,
... n_sp=0.6,
... threshold_sp=grid.node_x,
... discharge_field="surface_water__discharge")
>>> fr.run_one_step()
>>> sp.run_one_step(1.)
>>> z.reshape(grid.shape)[1, :] # doctest: +NORMALIZE_WHITESPACE
array([ 0. , 0.0647484 , 0.58634455, 2.67253503,
8.49212152, 20.92606987, 36. ])
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
Braun, J., Willett, S. (2013). A very efficient O(n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution. Geomorphology 180-181(C), 170-179. https://dx.doi.org/10.1016/j.geomorph.2012.10.008
Initialize the Fastscape stream power component. Note: a timestep, dt, can no longer be supplied to this component through the input file. It must instead be passed directly to the run method.
| Parameters: |
|
|---|
K¶Erodibility (units depend on m_sp).
__init__(grid, K_sp=0.001, m_sp=0.5, n_sp=1.0, threshold_sp=0.0, discharge_field='drainage_area', erode_flooded_nodes=True)[source]¶Initialize the Fastscape stream power component. Note: a timestep, dt, can no longer be supplied to this component through the input file. It must instead be passed directly to the run method.
| Parameters: |
|
|---|
run_one_step(dt)[source]¶Erode for a single time step.
This method implements the stream power erosion across one time interval, dt, following the Braun-Willett (2013) implicit Fastscape algorithm.
This follows Landlab standardized component design, and supercedes the
old driving method erode.
| Parameters: | dt (float) – Time-step size |
|---|
FireGenerator(grid, mean_fire_recurrence=1.0, shape_parameter=3.5, scale_parameter=None)[source]¶Bases: landlab.core.model_component.Component
Generate a random fire event or time series.
| Parameters: |
|
|---|
Generate a random fire event in time.
| Parameters: |
|
|---|
__init__(grid, mean_fire_recurrence=1.0, shape_parameter=3.5, scale_parameter=None)[source]¶Generate a random fire event in time.
| Parameters: |
|
|---|
generate_fire_recurrence()[source]¶Get time to next fire.
Finds the time to next fire (fire recurrence) based on the scale parameter (63.5% of fire Weibull distribution) and the shape parameter (describes the skew of the histogram, shape = 3.5 represents a normal distribution).
Rounds the time to next fire to 4 significant figures, for neatness.
| Returns: | Updated value for the time to next fire. |
|---|---|
| Return type: | float |
get_scale_parameter()[source]¶Get the scale parameter.
sets the scale parameter.
scale_parameter¶Scale parameter for the random distribution.
Flexure(grid, eet=65000.0, youngs=70000000000.0, method='airy', rho_mantle=3300.0, gravity=9.80665, n_procs=1)[source]¶Bases: landlab.core.model_component.Component
Deform the lithosphere with 1D or 2D flexure.
Landlab component that implements a 1 and 2D lithospheric flexure model.
Examples
>>> from landlab import RasterModelGrid
>>> from landlab.components.flexure import Flexure
>>> grid = RasterModelGrid((5, 4), xy_spacing=(1.e4, 1.e4))
>>> lith_press = grid.add_zeros(
... "lithosphere__overlying_pressure_increment", at="node"
... )
>>> flex = Flexure(grid)
>>> flex.name
'Flexure'
>>> flex.input_var_names
('lithosphere__overlying_pressure_increment',)
>>> flex.output_var_names
('lithosphere_surface__elevation_increment',)
>>> sorted(flex.units) # doctest: +NORMALIZE_WHITESPACE
[('lithosphere__overlying_pressure_increment', 'Pa'),
('lithosphere_surface__elevation_increment', 'm')]
>>> flex.grid.number_of_node_rows
5
>>> flex.grid.number_of_node_columns
4
>>> flex.grid is grid
True
>>> np.all(grid.at_node['lithosphere_surface__elevation_increment'] == 0.)
True
>>> np.all(grid.at_node['lithosphere__overlying_pressure_increment'] == 0.)
True
>>> flex.update()
>>> np.all(grid.at_node['lithosphere_surface__elevation_increment'] == 0.)
True
>>> load = grid.at_node['lithosphere__overlying_pressure_increment']
>>> load[4] = 1e9
>>> dz = grid.at_node['lithosphere_surface__elevation_increment']
>>> np.all(dz == 0.)
True
>>> flex.update()
>>> np.all(grid.at_node['lithosphere_surface__elevation_increment'] == 0.)
False
References
Required Software Citation(s) Specific to this Component
Hutton, E., Syvitski, J. (2008). Sedflux 2.0: An advanced process-response model that generates three-dimensional stratigraphy. Computers & Geosciences. 34(10), 1319-1337. https://dx.doi.org/10.1016/j.cageo.2008.02.013
Additional References
Lambeck, K.: Geophysical Geodesy, The Slow Deformations of the Earth, Clarendon Press, Oxford, UK, 718 pp., 1988.
Initialize the flexure component.
| Parameters: |
|
|---|
__init__(grid, eet=65000.0, youngs=70000000000.0, method='airy', rho_mantle=3300.0, gravity=9.80665, n_procs=1)[source]¶Initialize the flexure component.
| Parameters: |
|
|---|
alpha¶Flexure parameter (m).
eet¶Effective elastic thickness (m).
gamma_mantle¶Specific density of mantle (N/m^3).
gravity¶Acceleration due to gravity (m/s^2).
method¶Name of method used to calculate deflections.
rho_mantle¶Density of mantle (kg/m^3).
subside_loads(loads, out=None)[source]¶Subside surface due to multiple loads.
| Parameters: |
|
|---|---|
| Returns: | Deflections caused by the loading. |
| Return type: | ndarray of float |
youngs¶Young’s modulus of lithosphere (Pa).
Flexure1D(grid, eet=65000.0, youngs=70000000000.0, method='airy', rho_mantle=3300.0, rho_water=1030.0, gravity=9.80665, rows=None)[source]¶Bases: landlab.core.model_component.Component
Deform the lithosphere with 1D flexure.
Landlab component that implements a 1D lithospheric flexure model.
| Parameters: |
|
|---|
Examples
>>> from landlab import RasterModelGrid
>>> from landlab.components.flexure import Flexure1D
>>> grid = RasterModelGrid((5, 4), xy_spacing=(1.e4, 1.e4))
>>> lith_press = grid.add_zeros(
... "node",
... "lithosphere__increment_of_overlying_pressure")
>>> flex = Flexure1D(grid)
>>> flex.name
'1D Flexure Equation'
>>> flex.input_var_names
('lithosphere__increment_of_overlying_pressure',)
>>> flex.output_var_names
('lithosphere_surface__increment_of_elevation',)
>>> sorted(flex.units) # doctest: +NORMALIZE_WHITESPACE
[('lithosphere__increment_of_overlying_pressure', 'Pa'),
('lithosphere_surface__increment_of_elevation', 'm')]
>>> flex.grid.number_of_node_rows
5
>>> flex.grid.number_of_node_columns
4
>>> flex.grid is grid
True
>>> np.all(grid.at_node['lithosphere_surface__increment_of_elevation'] == 0.)
True
>>> np.all(grid.at_node['lithosphere__increment_of_overlying_pressure'] == 0.)
True
>>> flex.update()
>>> np.all(grid.at_node['lithosphere_surface__increment_of_elevation'] == 0.)
True
>>> load = grid.at_node['lithosphere__increment_of_overlying_pressure']
>>> load[4] = 1e9
>>> dz = grid.at_node['lithosphere_surface__increment_of_elevation']
>>> np.all(dz == 0.)
True
>>> flex.update()
>>> np.all(grid.at_node['lithosphere_surface__increment_of_elevation'] == 0.)
False
Initialize the flexure component.
| Parameters: |
|
|---|
__init__(grid, eet=65000.0, youngs=70000000000.0, method='airy', rho_mantle=3300.0, rho_water=1030.0, gravity=9.80665, rows=None)[source]¶Initialize the flexure component.
| Parameters: |
|
|---|
alpha¶Flexure parameter (m).
calc_flexure(x, loads, alpha, rigidity, out=None)[source]¶Subside surface due to multiple loads.
| Parameters: |
|
|---|---|
| Returns: | Deflections caused by the loading. |
| Return type: | ndarray of float |
dz_at_node¶eet¶Effective elastic thickness (m).
gamma_mantle¶Specific density of mantle (N/m^3).
gravity¶Acceleration due to gravity (m/s^2).
load_at_node¶method¶Name of method used to calculate deflections.
rho_mantle¶Density of mantle (kg/m^3).
rho_water¶Density of water (kg/m^3).
rigidity¶Flexural rigidity (N m).
subside_loads(loads, out=None)[source]¶Subside surface due to multiple loads.
| Parameters: |
|
|---|---|
| Returns: | Deflections caused by the loading. |
| Return type: | ndarray of float |
x_at_node¶youngs¶Young’s modulus of lithosphere (Pa).
FlowAccumulator(grid, surface='topographic__elevation', flow_director='FlowDirectorSteepest', runoff_rate=None, depression_finder=None, **kwargs)[source]¶Bases: landlab.core.model_component.Component
Component to accumulate flow and calculate drainage area.
This is accomplished by first finding flow directions by a user-specified method and then calculating the drainage area and discharge.
Optionally, spatially variable runoff can be set either by the model grid field ‘water__unit_flux_in’ or the input variable runoff_rate*.
Optionally a depression finding component can be specified and flow directing, depression finding, and flow routing can all be accomplished together.
NOTE: The perimeter nodes NEVER contribute to the accumulating flux, even if the gradients from them point inwards to the main body of the grid. This is because under Landlab definitions, perimeter nodes lack cells, so cannot accumulate any discharge.
FlowAccumulator stores as ModelGrid fields:
- Node array of drainage areas: ‘drainage_area’
- Node array of discharges: ‘surface_water__discharge’
- Node array containing downstream-to-upstream ordered list of node
- IDs: ‘flow__upstream_node_order’
- Node array of all but the first element of the delta data structure:
- flow__data_structure_delta. The first element is always zero.
The FlowDirector component will add additional ModelGrid fields. DirectToOne methods(Steepest/D4 and D8) and DirectToMany(DINF and MFD) use the same model grid field names. Some of these fields will be different shapes if a DirectToOne or a DirectToMany method is used.
The FlowDirectors store the following as ModelGrid fields:
- Node array of receivers (nodes that receive flow), or ITS OWN ID if
- there is no receiver: ‘flow__receiver_node’. This array is 2D for RouteToMany methods and has the shape (n-nodes x max number of receivers).
- Node array of flow proportions: ‘flow__receiver_proportions’. This
- array is 2D, for RouteToMany methods and has the shape (n-nodes x max number of receivers).
- Node array of links carrying flow: ‘flow__link_to_receiver_node’.
- This array is 2D for RouteToMany methods and has the shape (n-nodes x max number of receivers).
- Node array of downhill slopes from each receiver:
- ‘topographic__steepest_slope’ This array is 2D for RouteToMany methods and has the shape (n-nodes x max number of receivers).
- Boolean node array of all local lows: ‘flow__sink_flag’
- Link array identifing if flow goes with (1) or against (-1) the link direction: ‘flow__link_direction’
The primary method of this class is run_one_step.
run_one_step takes the optional argument update_flow_director (default is True) that determines if the flow_director is re-run before flow is accumulated.
| Parameters: |
|
|---|
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> mg = RasterModelGrid((3,3))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
The FlowAccumulator component accumulates flow and calculates drainage using all of the different methods for directing flow in Landlab. These include steepest descent (also known as D4 for the case of a raster grid) and D8 (on raster grids only). The method for flow director can be specified as a string (e.g., ‘D8’ or ‘FlowDirectorD8’), as an uninstantiated FlowDirector component or as an instantiated FlowDirector component.
The default method is to use FlowDirectorSteepest.
First let’s look at the three ways to instantiate a FlowAccumulator. The following four methods are all equivalent. First, we can pass the entire name of a flow director as a string to the argument flow_director:
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director='FlowDirectorSteepest'
... )
Second, we can pass just the method name as a string to the argument flow_director:
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director='Steepest'
... )
Third, we can import a FlowDirector component from Landlab and pass it to flow_director:
>>> from landlab.components import FlowDirectorSteepest
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
Finally, we can instantiate a FlowDirector component and pass this instantiated version to flow_director. You might want to do this if you used a FlowDirector in order to set up something before starting a time loop and then want to use the same flow director within the loop.
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
Now let’s look at what FlowAccumulator does. Even before we run FlowAccumulator it has the property surface_values that stores the values of the surface over which flow is directed and accumulated.
>>> fa.surface_values
array([ 0., 1., 2., 1., 2., 3., 2., 3., 4.])
Now let’s make a more complicated elevation grid for the next examples.
>>> mg = RasterModelGrid((5, 4))
>>> topographic__elevation = np.array([0., 0., 0., 0.,
... 0., 21., 10., 0.,
... 0., 31., 20., 0.,
... 0., 32., 30., 0.,
... 0., 0., 0., 0.])
>>> _ = mg.add_field("topographic__elevation", topographic__elevation, at="node")
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
>>> fa.run_one_step()
>>> mg.at_node['flow__receiver_node'] # doctest: +NORMALIZE_WHITESPACE
array([ 0, 1, 2, 3,
4, 1, 2, 7,
8, 10, 6, 11,
12, 14, 10, 15,
16, 17, 18, 19])
>>> mg.at_node['drainage_area'] # doctest: +NORMALIZE_WHITESPACE
array([ 0., 1., 5., 0.,
0., 1., 5., 0.,
0., 1., 4., 0.,
0., 1., 2., 0.,
0., 0., 0., 0.])
Now let’s change the cell area (100.) and the runoff rates:
>>> mg = RasterModelGrid((5, 4), xy_spacing=(10., 10))
Put the data back into the new grid.
>>> _ = mg.add_field("topographic__elevation", topographic__elevation, at="node")
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
>>> runoff_rate = np.arange(mg.number_of_nodes, dtype=float)
>>> rnff = mg.add_field("water__unit_flux_in", runoff_rate, at="node", clobber=True)
>>> fa.run_one_step()
>>> mg.at_node['surface_water__discharge'] # doctest: +NORMALIZE_WHITESPACE
array([ 0., 500., 5200., 0.,
0., 500., 5200., 0.,
0., 900., 4600., 0.,
0., 1300., 2700., 0.,
0., 0., 0., 0.])
The flow accumulator will happily work with a negative runoff rate, which could be used to allow, e.g., transmission losses:
>>> runoff_rate.fill(1.)
>>> fa.run_one_step()
>>> mg.at_node['surface_water__discharge']
array([ 0., 100., 500., 0.,
0., 100., 500., 0.,
0., 100., 400., 0.,
0., 100., 200., 0.,
0., 0., 0., 0.])
>>> runoff_rate[:8] = -0.5
>>> fa.run_one_step()
>>> mg.at_node['surface_water__discharge']
array([ 0., 0., 350., 0.,
0., 0., 350., 0.,
0., 100., 400., 0.,
0., 100., 200., 0.,
0., 0., 0., 0.])
The drainage area array is unaffected, as you would expect:
>>> mg.at_node['drainage_area']
array([ 0., 100., 500., 0.,
0., 100., 500., 0.,
0., 100., 400., 0.,
0., 100., 200., 0.,
0., 0., 0., 0.])
The FlowAccumulator component will work for both raster grids and irregular grids. For the example we will use a Hexagonal Model Grid, a special type of Voroni Grid that has regularly spaced hexagonal cells.
>>> from landlab import HexModelGrid
>>> hmg = HexModelGrid((5, 3), xy_of_lower_left=(-1., 0.))
>>> _ = hmg.add_field(
... "topographic__elevation",
... hmg.node_x + np.round(hmg.node_y),
... at="node",
... )
>>> fa = FlowAccumulator(
... hmg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
>>> fa.surface_values
array([ 0. , 1. , 2. ,
0.5, 1.5, 2.5, 3.5,
1. , 2. , 3. , 4. , 5. ,
2.5, 3.5, 4.5, 5.5,
3. , 4. , 5. ])
If the FlowDirector you want to use takes keyword arguments and you want to specify it using a string or uninstantiated FlowDirector class, include those keyword arguments when you create FlowAccumulator.
For example, in the case of a raster grid, FlowDirectorMFD can use only orthogonal links, or it can use both orthogonal and diagonal links.
>>> mg = RasterModelGrid((5, 5))
>>> topographic__elevation = mg.node_y+mg.node_x
>>> _ = mg.add_field("topographic__elevation", topographic__elevation, at="node")
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director='MFD',
... diagonals = True
... )
>>> fa.run_one_step()
>>> mg.at_node['flow__receiver_node'] # doctest: +NORMALIZE_WHITESPACE
array([[ 0, -1, -1, -1, -1, -1, -1, -1],
[ 1, -1, -1, -1, -1, -1, -1, -1],
[ 2, -1, -1, -1, -1, -1, -1, -1],
[ 3, -1, -1, -1, -1, -1, -1, -1],
[ 4, -1, -1, -1, -1, -1, -1, -1],
[ 5, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, 5, 1, -1, -1, 0, -1],
[-1, -1, 6, 2, -1, -1, 1, -1],
[-1, -1, 7, 3, -1, -1, 2, -1],
[ 9, -1, -1, -1, -1, -1, -1, -1],
[10, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, 10, 6, -1, -1, 5, -1],
[-1, -1, 11, 7, -1, -1, 6, -1],
[-1, -1, 12, 8, -1, -1, 7, -1],
[14, -1, -1, -1, -1, -1, -1, -1],
[15, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, 15, 11, -1, -1, 10, -1],
[-1, -1, 16, 12, -1, -1, 11, -1],
[-1, -1, 17, 13, -1, -1, 12, -1],
[19, -1, -1, -1, -1, -1, -1, -1],
[20, -1, -1, -1, -1, -1, -1, -1],
[21, -1, -1, -1, -1, -1, -1, -1],
[22, -1, -1, -1, -1, -1, -1, -1],
[23, -1, -1, -1, -1, -1, -1, -1],
[24, -1, -1, -1, -1, -1, -1, -1]])
>>> mg.at_node['drainage_area'].round(4) # doctest: +NORMALIZE_WHITESPACE
array([ 1.4117, 2.065 , 1.3254, 0.4038, 0. ,
2.065 , 3.4081, 2.5754, 1.3787, 0. ,
1.3254, 2.5754, 2.1716, 1.2929, 0. ,
0.4038, 1.3787, 1.2929, 1. , 0. ,
0. , 0. , 0. , 0. , 0. ])
It may seem odd that there are no round numbers in the drainage area field. This is because flow is directed to all downhill boundary nodes and partitioned based on slope.
To check that flow is conserved, sum along all boundary nodes.
>>> round(sum(mg.at_node['drainage_area'][mg.boundary_nodes]), 4)
9.0
This should be the same as the number of core nodes — as boundary nodes in landlab do not have area.
>>> len(mg.core_nodes)
9
Next, let’s set the dx spacing such that each cell has an area of one.
>>> dx=(2./(3.**0.5))**0.5
>>> hmg = HexModelGrid((5, 3), spacing=dx, xy_of_lower_left=(-1.0745, 0.))
>>> _ = hmg.add_field(
... "topographic__elevation",
... hmg.node_x**2 + np.round(hmg.node_y)**2,
... at="node",
... )
>>> fa = FlowAccumulator(
... hmg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
>>> fa.run_one_step()
>>> hmg.at_node['flow__receiver_node']
array([ 0, 1, 2,
3, 0, 1, 6,
7, 3, 4, 5, 11,
12, 8, 9, 15,
16, 17, 18])
>>> np.round(hmg.at_node['drainage_area'])
array([ 3., 2., 0.,
2., 3., 2., 0.,
0., 2., 2., 1., 0.,
0., 1., 1., 0.,
0., 0., 0.])
Now let’s change the cell area (100.) and the runoff rates:
>>> hmg = HexModelGrid((5, 3), spacing=dx * 10.0, xy_of_lower_left=(-10.745, 0.))
Put the data back into the new grid.
>>> _ = hmg.add_field(
... "topographic__elevation",
... hmg.node_x ** 2 + np.round(hmg.node_y) ** 2,
... at="node",
... )
>>> fa = FlowAccumulator(
... hmg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
>>> fa.run_one_step()
>>> np.round(hmg.at_node['surface_water__discharge'])
array([ 500., 0., 0.,
200., 500., 200., 0.,
0., 200., 200., 100., 0.,
0., 100., 100., 0.,
0., 0., 0.])
Next, let’s see what happens to a raster grid when there is a depression.
>>> mg = RasterModelGrid((7, 7), xy_spacing=0.5)
>>> z = mg.add_field("topographic__elevation", mg.node_x.copy(), at="node")
>>> z += 0.01 * mg.node_y
>>> mg.at_node['topographic__elevation'].reshape(mg.shape)[2:5, 2:5] *= 0.1
>>> mg.set_closed_boundaries_at_grid_edges(True, True, False, True)
This model grid has a depression in the center.
>>> mg.at_node['topographic__elevation'].reshape(mg.shape)
array([[ 0. , 0.5 , 1. , 1.5 , 2. , 2.5 , 3. ],
[ 0.005 , 0.505 , 1.005 , 1.505 , 2.005 , 2.505 , 3.005 ],
[ 0.01 , 0.51 , 0.101 , 0.151 , 0.201 , 2.51 , 3.01 ],
[ 0.015 , 0.515 , 0.1015, 0.1515, 0.2015, 2.515 , 3.015 ],
[ 0.02 , 0.52 , 0.102 , 0.152 , 0.202 , 2.52 , 3.02 ],
[ 0.025 , 0.525 , 1.025 , 1.525 , 2.025 , 2.525 , 3.025 ],
[ 0.03 , 0.53 , 1.03 , 1.53 , 2.03 , 2.53 , 3.03 ]])
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
>>> fa.run_one_step() # the flow "gets stuck" in the hole
>>> mg.at_node['flow__receiver_node'].reshape(mg.shape)
array([[ 0, 1, 2, 3, 4, 5, 6],
[ 7, 7, 16, 17, 18, 11, 13],
[14, 14, 16, 16, 17, 18, 20],
[21, 21, 16, 23, 24, 25, 27],
[28, 28, 23, 30, 31, 32, 34],
[35, 35, 30, 31, 32, 39, 41],
[42, 43, 44, 45, 46, 47, 48]])
>>> mg.at_node['drainage_area'].reshape(mg.shape)
array([[ 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
[ 0.25, 0.25, 0.25, 0.25, 0.5 , 0.25, 0. ],
[ 0.25, 0.25, 5. , 1.5 , 1. , 0.25, 0. ],
[ 0.25, 0.25, 3. , 0.75, 0.5 , 0.25, 0. ],
[ 0.25, 0.25, 2. , 1.5 , 1. , 0.25, 0. ],
[ 0.25, 0.25, 0.25, 0.25, 0.5 , 0.25, 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. , 0. ]])
Because of the depression, the flow ‘got stuck’ in the hole in the center of the grid. We can fix this by using a depression finder, such as DepressionFinderAndRouter.
>>> from landlab.components import DepressionFinderAndRouter
We can either run the depression finder separately from the flow accumulator or we can specify the depression finder and router when we instantiate the accumulator and it will run automatically. Similar to specifying the FlowDirector we can provide a depression finder in multiple three ways.
First let’s try running them separately.
>>> df_4 = DepressionFinderAndRouter(mg)
>>> df_4.map_depressions()
>>> mg.at_node['flow__receiver_node'].reshape(mg.shape)
array([[ 0, 1, 2, 3, 4, 5, 6],
[ 7, 7, 16, 17, 18, 11, 13],
[14, 14, 8, 16, 17, 18, 20],
[21, 21, 16, 16, 24, 25, 27],
[28, 28, 23, 24, 24, 32, 34],
[35, 35, 30, 31, 32, 39, 41],
[42, 43, 44, 45, 46, 47, 48]])
>>> mg.at_node['drainage_area'].reshape(mg.shape)
array([[ 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
[ 5.25, 5.25, 0.25, 0.25, 0.5 , 0.25, 0. ],
[ 0.25, 0.25, 5. , 1.5 , 1. , 0.25, 0. ],
[ 0.25, 0.25, 0.75, 2.25, 0.5 , 0.25, 0. ],
[ 0.25, 0.25, 0.5 , 0.5 , 1. , 0.25, 0. ],
[ 0.25, 0.25, 0.25, 0.25, 0.5 , 0.25, 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. , 0. ]])
Now the flow is routed correctly. The depression finder has properties that including whether there is a lake at the node, which lake is at each node, the outlet node of each lake, and the area of each lake.
>>> df_4.lake_at_node.reshape(mg.shape) # doctest: +NORMALIZE_WHITESPACE
array([[False, False, False, False, False, False, False],
[False, False, False, False, False, False, False],
[False, False, True, True, True, False, False],
[False, False, True, True, True, False, False],
[False, False, True, True, True, False, False],
[False, False, False, False, False, False, False],
[False, False, False, False, False, False, False]], dtype=bool)
>>> df_4.lake_map.reshape(mg.shape) # doctest: +NORMALIZE_WHITESPACE
array([[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, 16, 16, 16, -1, -1],
[-1, -1, 16, 16, 16, -1, -1],
[-1, -1, 16, 16, 16, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1]])
>>> df_4.lake_codes # a unique code for each lake present on the grid
array([16])
>>> df_4.lake_outlets # the outlet node of each lake in lake_codes
array([8])
>>> df_4.lake_areas # the area of each lake in lake_codes
array([ 2.25])
Alternatively, we can initialize a flow accumulator with a depression finder specified. Calling run_one_step() will run both the accumulator and the depression finder with one call. For this example, we will pass the class DepressionFinderAndRouter to the parameter depression_finder.
>>> mg = RasterModelGrid((7, 7), xy_spacing=0.5)
>>> z = mg.add_field("topographic__elevation", mg.node_x.copy(), at="node")
>>> z += 0.01 * mg.node_y
>>> mg.at_node['topographic__elevation'].reshape(mg.shape)[2:5, 2:5] *= 0.1
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director='FlowDirectorD8',
... depression_finder=DepressionFinderAndRouter
... )
>>> fa.run_one_step()
This has the same effect of first calling the accumulator and then calling the depression finder.
>>> mg.at_node['flow__receiver_node'].reshape(mg.shape)
array([[ 0, 1, 2, 3, 4, 5, 6],
[ 7, 7, 16, 17, 18, 18, 13],
[14, 14, 8, 16, 17, 18, 20],
[21, 21, 16, 16, 24, 25, 27],
[28, 28, 23, 24, 24, 32, 34],
[35, 35, 30, 31, 32, 32, 41],
[42, 43, 44, 45, 46, 47, 48]])
>>> mg.at_node['drainage_area'].reshape(mg.shape)
array([[ 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
[ 5.25, 5.25, 0.25, 0.25, 0.25, 0.25, 0. ],
[ 0.25, 0.25, 5. , 1.5 , 1. , 0.25, 0. ],
[ 0.25, 0.25, 0.75, 2.25, 0.5 , 0.25, 0. ],
[ 0.25, 0.25, 0.5 , 0.5 , 1. , 0.25, 0. ],
[ 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. , 0. ]])
The depression finder is stored as part of the flow accumulator, so its properties can be accessed through the depression finder.
>>> fa.depression_finder.lake_at_node.reshape(mg.shape) # doctest: +NORMALIZE_WHITESPACE
array([[False, False, False, False, False, False, False],
[False, False, False, False, False, False, False],
[False, False, True, True, True, False, False],
[False, False, True, True, True, False, False],
[False, False, True, True, True, False, False],
[False, False, False, False, False, False, False],
[False, False, False, False, False, False, False]], dtype=bool)
>>> fa.depression_finder.lake_map.reshape(mg.shape) # doctest: +NORMALIZE_WHITESPACE
array([[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, 16, 16, 16, -1, -1],
[-1, -1, 16, 16, 16, -1, -1],
[-1, -1, 16, 16, 16, -1, -1],
[-1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1]])
>>> fa.depression_finder.lake_codes # a unique code for each lake present on the grid
array([16])
>>> fa.depression_finder.lake_outlets # the outlet node of each lake in lake_codes
array([8])
>>> fa.depression_finder.lake_areas # the area of each lake in lake_codes
array([ 2.25])
Finally, note that the DepressionFinderAndRouter takes a keyword argument routing (‘D8’, default; ‘D4’) that sets how connectivity is set between nodes. Similar to our ability to pass keyword arguments to the FlowDirector through FlowAccumulator, we can pass this keyword argument to the DepressionFinderAndRouter component.
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest,
... depression_finder=DepressionFinderAndRouter,
... routing='D4'
... )
FlowAccumulator was designed to work with all types of grids. However,
NetworkModelGrid’s have no cell area. Thus, in order for FlowAccumulator to
this type of grid, an at-node array called cell_area_at_node must be
present.
>>> from landlab.grid.network import NetworkModelGrid
>>> y_of_node = (0, 1, 2, 2)
>>> x_of_node = (0, 0, -1, 1)
>>> nodes_at_link = ((1, 0), (2, 1), (3, 1))
>>> nmg = NetworkModelGrid((y_of_node, x_of_node), nodes_at_link)
>>> area = nmg.add_ones("cell_area_at_node", at="node")
>>> z = nmg.add_field(
... "topographic__elevation",
... nmg.x_of_node + nmg.y_of_node,
... at="node",
... )
>>> fa = FlowAccumulator(nmg)
>>> fa.run_one_step()
>>> nmg.at_node['flow__receiver_node']
array([0, 0, 2, 1])
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
Braun, J., Willett, S. (2013). A very efficient O(n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution. Geomorphology 180-181(C), 170-179. https://dx.doi.org/10.1016/j.geomorph.2012.10.008
Initialize the FlowAccumulator component.
Saves the grid, tests grid type, tests imput types and compatability for the flow_director and depression_finder keyword arguments, tests the argument of runoff_rate, and initializes new fields.
__init__(grid, surface='topographic__elevation', flow_director='FlowDirectorSteepest', runoff_rate=None, depression_finder=None, **kwargs)[source]¶Initialize the FlowAccumulator component.
Saves the grid, tests grid type, tests imput types and compatability for the flow_director and depression_finder keyword arguments, tests the argument of runoff_rate, and initializes new fields.
accumulate_flow(update_flow_director=True, update_depression_finder=True)[source]¶Function to make FlowAccumulator calculate drainage area and discharge.
Running accumulate_flow() results in the following to occur:
- Flow directions are updated (unless update_flow_director is set as False). This incluldes checking for updated boundary conditions.
- The depression finder, if present is updated (unless update_depression_finder is set as False).
- Intermediate steps that analyse the drainage network topology and create datastructures for efficient drainage area and discharge calculations.
- Calculation of drainage area and discharge.
- Return of drainage area and discharge.
| Parameters: |
|
|---|---|
| Returns: |
|
depression_finder¶The DepressionFinder used internally.
flow_director¶The FlowDirector used internally.
headwater_nodes()[source]¶Return the headwater nodes.
These are nodes that contribute flow and have no upstream nodes.
Examples
>>> from numpy.testing import assert_array_equal
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> mg = RasterModelGrid((5, 5))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
>>> fa = FlowAccumulator(mg, 'topographic__elevation')
>>> fa.run_one_step()
>>> assert_array_equal(fa.headwater_nodes(), np.array([16, 17, 18]))
link_order_upstream()[source]¶Return the upstream order of active links.
Examples
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> mg = RasterModelGrid((5, 5))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
>>> fa = FlowAccumulator(mg, 'topographic__elevation')
>>> fa.run_one_step()
>>> fa.link_order_upstream()
array([ 5, 14, 23, 6, 15, 24, 7, 16, 25])
This also works for route-to-many methods
>>> mg = RasterModelGrid((5, 5))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
>>> fa = FlowAccumulator(mg,
... 'topographic__elevation',
... flow_director='MFD')
>>> fa.run_one_step()
>>> fa.link_order_upstream()
array([ 5, 14, 10, 6, 11, 7, 23, 19, 15, 20, 16, 28, 24, 29, 25])
node_drainage_area¶Return the drainage area.
node_order_upstream¶Return the upstream node order (drainage stack).
node_water_discharge¶Return the surface water discharge.
run_one_step()[source]¶Accumulate flow and save to the model grid.
An alternative to run_one_step() is accumulate_flow() which does the same things but also returns the drainage area and discharge. accumulate_flow() additionally provides the ability to turn off updating the flow director or the depression finder.
surface_values¶Values of the surface over which flow is accumulated.
FlowDirectorD8(grid, surface='topographic__elevation')[source]¶Bases: landlab.components.flow_director.flow_director_to_one._FlowDirectorToOne
Single-path (steepest direction) flow direction with diagonals on rasters.
Single-path (steepest direction) flow direction finding on raster grids by the D8 method. This method considers flow on all eight links such that flow is possible on orthogonal and on diagonal links.
The method that considers only orthogonal links (D4 method) for raster grids is FlowDirectorSteepest.
This method is not implemented for Voroni grids, use FlowDirectorSteepest instead.
Stores as ModelGrid fields:
The primary method of this class is run_one_step.
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorD8
>>> mg = RasterModelGrid((3,3), xy_spacing=(1, 1))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
>>> fd = FlowDirectorD8(mg, 'topographic__elevation')
>>> fd.surface_values
array([ 0., 1., 2., 1., 2., 3., 2., 3., 4.])
>>> fd.run_one_step()
>>> mg.at_node['flow__receiver_node']
array([0, 1, 2, 3, 0, 5, 6, 7, 8])
>>> mg.at_node['topographic__steepest_slope']
array([ 0. , 0. , 0. , 0. , 1.41421356,
0. , 0. , 0. , 0. ])
>>> mg.at_node['flow__link_to_receiver_node']
array([-1, -1, -1, -1, 12, -1, -1, -1, -1])
>>> mg.at_node['flow__sink_flag'].astype(int)
array([1, 1, 1, 1, 0, 1, 1, 1, 1])
>>> mg_2 = RasterModelGrid((5, 4), xy_spacing=(1, 1))
>>> topographic__elevation = np.array([0., 0., 0., 0.,
... 0., 21., 10., 0.,
... 0., 31., 20., 0.,
... 0., 32., 30., 0.,
... 0., 0., 0., 0.])
>>> _ = mg_2.add_field(
... "topographic__elevation",
... topographic__elevation,
... at="node",
... )
>>> mg_2.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fd_2 = FlowDirectorD8(mg_2)
>>> fd_2.run_one_step()
>>> mg_2.at_node['flow__receiver_node'] # doctest: +NORMALIZE_WHITESPACE
array([ 0, 1, 2, 3,
4, 1, 2, 7,
8, 6, 6, 11,
12, 10, 10, 15,
16, 17, 18, 19])
The flow directors also have the ability to return the flow receiver nodes
>>> receiver = fd.direct_flow()
>>> receiver
array([0, 1, 2,
3, 0, 5,
6, 7, 8])
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
O’Callaghan, J., Mark, D. (1984). The extraction of drainage networks from digital elevation data. Computer Vision, Graphics, and Image Processing 28(3), 323 - 344. https://dx.doi.org/10.1016/s0734-189x(84)80011-0
| Parameters: |
|
|---|
__init__(grid, surface='topographic__elevation')[source]¶| Parameters: |
|
|---|
direct_flow()[source]¶Find flow directions, save to the model grid, and return receivers.
direct_flow() checks for updated boundary conditions, calculates slopes on links, finds baselevel nodes based on the status at node, calculates flow directions, saves results to the grid, and returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_node’]
an alternative to direct_flow() is run_one_step() which does the same things but also returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_node’]
run_one_step()[source]¶Find flow directions and save to the model grid.
run_one_step() checks for updated boundary conditions, calculates slopes on links, finds baselevel nodes based on the status at node, calculates flow directions, and saves results to the grid.
an alternative to direct_flow() is direct_flow() which does the same things but also returns the receiver nodes not return values.
FlowDirectorDINF(grid, surface='topographic__elevation')[source]¶Bases: landlab.components.flow_director.flow_director_to_many._FlowDirectorToMany
Flow direction on a raster grid by the D infinity method.
Directs flow by the D infinity method (Tarboton, 1997). Each node is assigned two flow directions, toward the two neighboring nodes that are on the steepest subtriangle. Partitioning of flow is done based on the aspect of the subtriangle.
Specifically, it stores as ModelGrid fields:
The primary method of this class is run_one_step.
Examples
This method works for both raster and irregular grids. First we will look at a raster example, and then an irregular example.
>>> import numpy as numpy
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorDINF
>>> mg = RasterModelGrid((4,4), xy_spacing=(1, 1))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x**2 + mg.node_y**2,
... at="node",
... )
The DINF flow director can be uses for raster grids only.
>>> fd = FlowDirectorDINF(mg, 'topographic__elevation')
>>> fd.surface_values.reshape(mg.shape)
array([[ 0., 1., 4., 9.],
[ 1., 2., 5., 10.],
[ 4., 5., 8., 13.],
[ 9., 10., 13., 18.]])
>>> fd.run_one_step()
Unlike flow directors that only direct flow to one node or to all downstream nodes, FlowDirectorDINF directs flow two nodes only. It stores the receiver information is a (number of nodes x 2) shape field at nodes.
>>> mg.at_node['flow__receiver_node']
array([[ 0, -1],
[ 1, -1],
[ 2, -1],
[ 3, -1],
[ 4, -1],
[ 0, -1],
[ 5, 1],
[ 7, -1],
[ 8, -1],
[ 5, -1],
[ 5, -1],
[11, -1],
[12, -1],
[13, -1],
[14, -1],
[15, -1]])
It also stores the proportions of flow going to each receiver, the link on which the flow moves in at node arrays, and the slope of each link.
>>> mg.at_node['flow__receiver_proportions'] # doctest: +NORMALIZE_WHITESPACE
array([[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 0.59033447, 0.40966553],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ]])
>>> mg.at_node['flow__link_to_receiver_node']
array([[-1, -1],
[-1, -1],
[-1, -1],
[-1, -1],
[ 3, 25],
[24, 4],
[ 8, 26],
[ 9, 28],
[14, 31],
[11, 33],
[32, 12],
[16, 34],
[21, 37],
[18, 39],
[19, 38],
[20, 40]])
>>> mg.at_node['topographic__steepest_slope'] # doctest: +NORMALIZE_WHITESPACE
array([[ -1.00000000e+00, -1.41421356e+00],
[ 1.00000000e+00, -7.12763635e+02],
[ 3.00000000e+00, 1.41421356e+00],
[ 5.00000000e+00, 2.82842712e+00],
[ 1.00900000e+03, 7.12763635e+02],
[ 1.41421356e+00, 1.00000000e+00],
[ 3.00000000e+00, 2.82842712e+00],
[ 1.00400000e+03, 7.10642315e+02],
[ 1.00400000e+03, 7.12056529e+02],
[ 3.00000000e+00, 0.00000000e+00],
[ 4.24264069e+00, 3.00000000e+00],
[ 1.00100000e+03, 7.09935208e+02],
[ -0.00000000e+00, 7.09935208e+02],
[ 1.00400000e+03, 7.07813888e+02],
[ 1.00100000e+03, 7.09935208e+02],
[ 0.00000000e+00, 7.07813888e+02]])
Finally, FlowDirectorDINF identifies sinks, or local lows.
>>> mg.at_node['flow__sink_flag'].astype(int)
array([1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1])
The flow directors also have the ability to return the flow receiver nodes through a function called direct_flow()
>>> fd = FlowDirectorDINF(mg, 'topographic__elevation')
>>> fd.run_one_step()
>>> receivers, proportions = fd.direct_flow()
>>> receivers
array([[ 0, -1],
[ 1, -1],
[ 2, -1],
[ 3, -1],
[ 4, -1],
[ 0, -1],
[ 5, 1],
[ 7, -1],
[ 8, -1],
[ 5, -1],
[ 5, -1],
[11, -1],
[12, -1],
[13, -1],
[14, -1],
[15, -1]])
>>> proportions # doctest: +NORMALIZE_WHITESPACE
array([[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 0.59033447, 0.40966553],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ],
[ 1. , 0. ]])
For each donor node (represented by each row) the proportions should sum to one.
>>> proportions.sum(axis=1)
array([ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1.])
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
Tarboton, D. (1997). A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resources Research 33(2), 309-319. https://dx.doi.org/10.1029/96wr03137
| Parameters: |
|
|---|
__init__(grid, surface='topographic__elevation')[source]¶| Parameters: |
|
|---|
direct_flow()[source]¶Find flow directions, save to the model grid, and return receivers.
direct_flow() checks for updated boundary conditions, calculates slopes on links, finds basself._surface_valuesel nodes based on the status at node, calculates flow directions, saves results to the grid, and returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_nodes’]
An alternative to direct_flow() is run_one_step() which does the same things but also returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_nodes’]
run_one_step()[source]¶Find flow directions and save to the model grid.
run_one_step() checks for updated boundary conditions, calculates slopes on links, finds basself._surface_valuesel nodes based on the status at node, calculates flow directions, and saves results to the grid.
An alternative to direct_flow() is direct_flow() which does the same things but also returns the receiver nodes not return values.
FlowDirectorMFD(grid, surface='topographic__elevation', **kwargs)[source]¶Bases: landlab.components.flow_director.flow_director_to_many._FlowDirectorToMany
Multiple-path flow direction with or without out diagonals.
Directs flow by the multiple flow direction method. Each node is assigned multiple flow directions, toward all of the N neighboring nodes that are lower than it. If none of the neighboring nodes are lower, the location is identified as a pit. Flow proportions can be calculated as proportional to slope or proportional to the square root of slope, which is the solution to a steady kinematic wave.
Specifically, it stores as ModelGrid fields:
The primary method of this class is run_one_step.
Examples
This method works for both raster and irregular grids. First we will look at a raster example, and then an irregular example.
>>> import numpy as numpy
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorMFD
>>> mg = RasterModelGrid((3,3), xy_spacing=(1, 1))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
The MFD flow director can be uses for raster and irregular grids. For raster grids, use of diagonal links is specified with the keyword diagonals (default is False).
>>> fd = FlowDirectorMFD(mg, 'topographic__elevation', diagonals = True)
>>> fd.surface_values
array([ 0., 1., 2., 1., 2., 3., 2., 3., 4.])
>>> fd.run_one_step()
Unlike flow directors that only direct flow to one node, FlowDirectorMFD directs flow to all downstream nodes. It stores the receiver information is a (number of nodes x maximum number or receivers) shape field at nodes.
>>> mg.at_node['flow__receiver_node']
array([[ 0, -1, -1, -1, -1, -1, -1, -1],
[ 1, -1, -1, -1, -1, -1, -1, -1],
[ 2, -1, -1, -1, -1, -1, -1, -1],
[ 3, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, 1, -1, -1, 0, -1],
[ 5, -1, -1, -1, -1, -1, -1, -1],
[ 6, -1, -1, -1, -1, -1, -1, -1],
[ 7, -1, -1, -1, -1, -1, -1, -1],
[ 8, -1, -1, -1, -1, -1, -1, -1]])
It also stores the proportions of flow going to each receiver, the link on which the flow moves in at node arrays, and the slope of each link.
>>> mg.at_node['flow__receiver_proportions'] # doctest: +NORMALIZE_WHITESPACE
array([[ 1. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ],
[ 0. , 0. , 0. , 0.41421356, 0. ,
0. , 0.58578644, 0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ]])
>>> mg.at_node['flow__link_to_receiver_node']
array([[-1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, 3, -1, -1, 12, -1],
[-1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1, -1]])
>>> mg.at_node['topographic__steepest_slope'] # doctest: +NORMALIZE_WHITESPACE
array([[ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ],
[ 0. , 0. , 0. , 1. , 0. ,
0. , 1.41421356, 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ]])
Finally, FlowDirectorMFD identifies sinks, or local lows.
>>> mg.at_node['flow__sink_flag'].astype(int)
array([1, 1, 1, 1, 0, 1, 1, 1, 1])
The flow directors also have the ability to return the flow receiver nodes. For this example, we will turn the diagonals off. This is the default value.
>>> mg = RasterModelGrid((3,3), xy_spacing=(1, 1))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
>>> fd = FlowDirectorMFD(mg, 'topographic__elevation')
>>> fd.run_one_step()
>>> receivers, proportions = fd.direct_flow()
>>> receivers
array([[ 0, -1, -1, -1],
[ 1, -1, -1, -1],
[ 2, -1, -1, -1],
[ 3, -1, -1, -1],
[-1, -1, -1, 1],
[ 5, -1, -1, -1],
[ 6, -1, -1, -1],
[ 7, -1, -1, -1],
[ 8, -1, -1, -1]])
>>> proportions # doctest: +NORMALIZE_WHITESPACE
array([[ 1., 0., 0., 0.],
[ 1., 0., 0., 0.],
[ 1., 0., 0., 0.],
[ 1., 0., 0., 0.],
[ 0., 0., 0., 1.],
[ 1., 0., 0., 0.],
[ 1., 0., 0., 0.],
[ 1., 0., 0., 0.],
[ 1., 0., 0., 0.]])
For each donor node (represented by each row) the proportions should sum to one.
>>> proportions.sum(axis=1)
array([ 1., 1., 1., 1., 1., 1., 1., 1., 1.])
For the second example we will use a Hexagonal Model Grid, a special type of Voroni Grid that has regularly spaced hexagonal cells. FlowDirectorMFD has multiple ways to partition flow based on slope. The default method is based on the slope angle. A secondary methods is to partion based on the square root of slope. This represents the solution to a steady kinematic wave.
>>> from landlab import HexModelGrid
>>> mg = HexModelGrid((5, 3))
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + numpy.round(mg.node_y),
... at="node",
... )
>>> fd = FlowDirectorMFD(
... mg,
... 'topographic__elevation',
... partition_method='square_root_of_slope'
... )
>>> fd.surface_values # doctest: +NORMALIZE_WHITESPACE
array([ 1. , 2. , 3. ,
1.5, 2.5, 3.5, 4.5,
2. , 3. , 4. , 5. , 6. ,
3.5, 4.5, 5.5, 6.5,
4. , 5. , 6. ])
>>> fd.run_one_step()
>>> mg.at_node['flow__receiver_node']
array([[ 0, -1, -1, -1, -1, -1],
[ 1, -1, -1, -1, -1, -1],
[ 2, -1, -1, -1, -1, -1],
[ 3, -1, -1, -1, -1, -1],
[-1, -1, -1, 3, 0, 1],
[-1, -1, -1, 4, 1, 2],
[ 6, -1, -1, -1, -1, -1],
[ 7, -1, -1, -1, -1, -1],
[-1, -1, -1, 7, 3, 4],
[-1, -1, -1, 8, 4, 5],
[-1, -1, -1, 9, 5, 6],
[11, -1, -1, -1, -1, -1],
[12, -1, -1, -1, -1, -1],
[-1, -1, 16, 12, 8, 9],
[-1, -1, 17, 13, 9, 10],
[15, -1, -1, -1, -1, -1],
[16, -1, -1, -1, -1, -1],
[17, -1, -1, -1, -1, -1],
[18, -1, -1, -1, -1, -1]])
>>> mg.at_node['flow__receiver_proportions'] # doctest: +NORMALIZE_WHITESPACE
array([[ 1. , 0. , 0. , 0. , 0. ,
0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. ],
[ 0. , 0. , 0. , 0.34108138, 0.41773767,
0.24118095],
[ 0. , 0. , 0. , 0.34108138, 0.41773767,
0.24118095],
[ 1. , 0. , 0. , 0. , 0. ,
0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. ],
[ 0. , 0. , 0. , 0.34108138, 0.41773767,
0.24118095],
[ 0. , 0. , 0. , 0.34108138, 0.41773767,
0.24118095],
[ 0. , 0. , 0. , 0.34108138, 0.41773767,
0.24118095],
[ 1. , 0. , 0. , 0. , 0. ,
0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. ],
[ 0. , 0. , 0.19431571, 0.27480391, 0.33656468,
0.19431571],
[ 0. , 0. , 0.19431571, 0.27480391, 0.33656468,
0.19431571],
[ 1. , 0. , 0. , 0. , 0. ,
0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. ],
[ 1. , 0. , 0. , 0. , 0. ,
0. ]])
>>> mg.at_node['flow__link_to_receiver_node']
array([[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, 8, 3, 4],
[-1, -1, -1, 9, 5, 6],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, 19, 12, 13],
[-1, -1, -1, 20, 14, 15],
[-1, -1, -1, 21, 16, 17],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1],
[-1, -1, 35, 31, 25, 26],
[-1, -1, 37, 32, 27, 28],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1]])
>>> mg.at_node['topographic__steepest_slope'] # doctest: +NORMALIZE_WHITESPACE
array([[ 0. , 0. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 1. , 1.5, 0.5],
[ 0. , 0. , 0. , 1. , 1.5, 0.5],
[ 0. , 0. , 1. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 1. , 1.5, 0.5],
[ 0. , 0. , 0. , 1. , 1.5, 0.5],
[ 0. , 0. , 0. , 1. , 1.5, 0.5],
[ 0. , 1. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0.5, 0. , 0. ],
[ 0. , 0. , 0.5, 1. , 1.5, 0.5],
[ 0. , 0. , 0.5, 1. , 1.5, 0.5],
[ 0. , 1. , 1.5, 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0.5, 0. , 0. , 0. ],
[ 0. , 0.5, 0. , 0. , 0. , 0. ]])
>>> mg.at_node['flow__sink_flag'].astype(int)
array([1, 1, 1,
1, 0, 0, 1,
1, 0, 0, 0, 1,
1, 0, 0, 1,
1, 1, 1])
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
Freeman, T. (1991). Calculating catchment area with divergent flow based on a regular grid. Computers and Geosciences 17(3), 413 - 422. https://dx.doi.org/10.1016/0098-3004(91)90048-i
Quinn, P., Beven, K., Chevallier, P., Planchon, O. (1991). The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models Hydrological Processes 5(1), 59-79. https://dx.doi.org/10.1002/hyp.3360050106
| Parameters: |
|
|---|
__init__(grid, surface='topographic__elevation', **kwargs)[source]¶| Parameters: |
|
|---|
direct_flow()[source]¶Find flow directions, save to the model grid, and return receivers.
direct_flow() checks for updated boundary conditions, calculates slopes on links, finds basself._surface_valuesel nodes based on the status at node, calculates flow directions, saves results to the grid, and returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_nodes’]
An alternative to direct_flow() is run_one_step() which does the same things but also returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_nodes’]
run_one_step()[source]¶Find flow directions and save to the model grid.
run_one_step() checks for updated boundary conditions, calculates slopes on links, finds basself._surface_valuesel nodes based on the status at node, calculates flow directions, and saves results to the grid.
An alternative to run_one_step() is direct_flow() which does the same things but also returns the receiver nodes not return values.
FlowDirectorSteepest(grid, surface='topographic__elevation')[source]¶Bases: landlab.components.flow_director.flow_director_to_one._FlowDirectorToOne
Single-path (steepest direction) flow direction without diagonals.
This components finds the steepest single-path steepest descent flow directions. It is equivalent to D4 method in the special case of a raster grid in that it does not consider diagonal links between nodes. For that capability, use FlowDirectorD8.
Stores as ModelGrid fields:
The primary method of this class is run_one_step.
Examples
This method works for both raster and irregular grids. First we will look at a raster example, and then an irregular example.
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorSteepest
>>> mg = RasterModelGrid((3,3), xy_spacing=(1, 1))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fd.surface_values
array([ 0., 1., 2., 1., 2., 3., 2., 3., 4.])
>>> fd.run_one_step()
>>> mg.at_node['flow__receiver_node']
array([0, 1, 2, 3, 1, 5, 6, 7, 8])
>>> mg.at_node['topographic__steepest_slope']
array([ 0., 0., 0., 0., 1., 0., 0., 0., 0.])
>>> mg.at_node['flow__link_to_receiver_node']
array([-1, -1, -1, -1, 3, -1, -1, -1, -1])
>>> mg.at_node['flow__sink_flag'].astype(int)
array([1, 1, 1, 1, 0, 1, 1, 1, 1])
>>> mg_2 = RasterModelGrid((5, 4), xy_spacing=(1, 1))
>>> topographic__elevation = np.array([0., 0., 0., 0.,
... 0., 21., 10., 0.,
... 0., 31., 20., 0.,
... 0., 32., 30., 0.,
... 0., 0., 0., 0.])
>>> _ = mg_2.add_field(
... "topographic__elevation",
... topographic__elevation,
... at="node",
... )
>>> mg_2.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fd_2 = FlowDirectorSteepest(mg_2)
>>> fd_2.run_one_step()
>>> mg_2.at_node['flow__receiver_node'] # doctest: +NORMALIZE_WHITESPACE
array([ 0, 1, 2, 3,
4, 1, 2, 7,
8, 10, 6, 11,
12, 14, 10, 15,
16, 17, 18, 19])
And the at-link field 'flow__link_direction' indicates if the flow along
the link is with or against the direction indicated by 'link_dirs_at_node'
(from tail node to head node).
>>> mg_2.at_link['flow__link_direction']
array([ 0, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0,
0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], dtype=int8)
This indicates that flow on links 4, 5, 12, and 19 goes against the topologic ordering – that is that flow goes from head node to tail node – and that flow goes with the topologic ordering on links 15 and 22. All other links have no flow on them.
The FlowDirectorSteepest attribute flow_link_direction_at_node indicates
the link flow direction (with or against topology directions) for all links
at node. The ordering of links at node mirrors the grid attribute
links_at_node.
>>> fd_2.flow_link_direction_at_node()
array([[ 0, 0, 0, 0],
[ 0, -1, 0, 0],
[ 0, -1, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, -1],
[ 0, -1, 0, -1],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 1, 0, 0, 0],
[ 0, -1, 1, -1],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 1, 0, 0, 0],
[ 0, 0, 1, -1],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0]], dtype=int8)
For example, this indicates that node 10 has flow going along three links that are attached to it. The link to the East has no flow, the link to the North has flow going against the topologic direction, the link to the West has flow going with the topologic direction, and the link to the South has flow going against the topologic direction.
In many use cases, one might want to know which links are bringing flow into
or out of the node. The flow director attribute flow_link_incoming_at_node
provides this information. Here -1 means that flow is outgoing from the node
and 1 means it is incoming.
>>> fd_2.flow_link_incoming_at_node()
array([[ 0, 0, 0, 0],
[ 0, 1, 0, 0],
[ 0, 1, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, -1],
[ 0, 1, 0, -1],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[-1, 0, 0, 0],
[ 0, 1, 1, -1],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[-1, 0, 0, 0],
[ 0, 0, 1, -1],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0]], dtype=int8)
So if one wanted to identify the source nodes at node, you would do the following:
>>> np.where(fd_2.flow_link_incoming_at_node() == 1, mg_2.adjacent_nodes_at_node, -1)
array([[-1, -1, -1, -1],
[-1, 5, -1, -1],
[-1, 6, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, 10, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, 14, 9, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, -1, 13, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1]])
The flow directors also have the ability to return the flow receiver nodes
>>> receiver = fd.direct_flow()
>>> receiver
array([0, 1, 2,
3, 1, 5,
6, 7, 8])
For the second example we will use a Hexagonal Model Grid, a special type of Voroni Grid that has regularly spaced hexagonal cells.
>>> from landlab import HexModelGrid
>>> mg = HexModelGrid((5, 3))
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + np.round(mg.node_y),
... at="node",
... )
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fd.surface_values
array([ 1. , 2. , 3. ,
1.5, 2.5, 3.5, 4.5,
2. , 3. , 4. , 5. , 6. ,
3.5, 4.5, 5.5, 6.5,
4. , 5. , 6. ])
>>> fd.run_one_step()
>>> mg.at_node['flow__receiver_node']
array([ 0, 1, 2,
3, 0, 1, 6,
7, 3, 4, 5, 11,
12, 8, 9, 15,
16, 17, 18])
>>> mg.at_node['topographic__steepest_slope']
array([ 0. , 0. , 0. ,
0. , 1.5, 1.5, 0. ,
0. , 1.5, 1.5, 1.5, 0. ,
0. , 1.5, 1.5, 0. ,
0. , 0. , 0. ])
>>> mg.at_node['flow__link_to_receiver_node']
array([-1, -1, -1,
-1, 3, 5, -1,
-1, 12, 14, 16, -1,
-1, 25, 27, -1,
-1, -1, -1])
>>> mg.at_node['flow__sink_flag'].astype(int)
array([1, 1, 1,
1, 0, 0, 1,
1, 0, 0, 0, 1,
1, 0, 0, 1,
1, 1, 1])
>>> receiver = fd.direct_flow()
>>> receiver
array([ 0, 1, 2,
3, 0, 1, 6,
7, 3, 4, 5, 11,
12, 8, 9, 15,
16, 17, 18])
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
None Listed
| Parameters: |
|
|---|
__init__(grid, surface='topographic__elevation')[source]¶| Parameters: |
|
|---|
direct_flow()[source]¶Find flow directions, save to the model grid, and return receivers.
direct_flow() checks for updated boundary conditions, calculates slopes on links, finds baselevel nodes based on the status at node, calculates flow directions, saves results to the grid, and returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_node’]
An alternative to direct_flow() is run_one_step() which does the same things but also returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_node’]
downstream_node_at_link()[source]¶At-link array of the downstream node based on flow direction.
BAD_INDEX_VALUE is given if no downstream node is defined.
Examples
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorSteepest
>>> mg = RasterModelGrid((3,3))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fd.run_one_step()
>>> fd.downstream_node_at_link()
array([-1, -1, -1, 1, -1, -1, -1, -1, -1, -1, -1, -1])
flow_link_direction¶Return array of flow link direction.
This property indicates the relationship between the flow direction (determined based on the elevation of nodes) and the topologic link direction (in which the head and tail nodes are defined based on relative position in x-y space).
It has the shape (number_of_links,).
Recall that the standard landlab link direction goes from the tail node to the head node.
A value of zero indicates that the link does not exist or is not active.
A value of -1 indicates that water flow based on
flow__link_to_receiver_node goes from head node to tail node, while
a value of 1 indicates that water flow goes from tail node to head
node.
Examples
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorSteepest
>>> mg = RasterModelGrid((3,3))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fd.run_one_step()
>>> fd.flow_link_direction
array([ 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0], dtype=int8)
flow_link_direction_at_node()[source]¶Return array of flow link direction at node.
This property mirrors links_at_node and indicates the relationship between the flow direction (determined based on the elevation of nodes) and the topologic link direction (in which the head and tail nodes are defined based on relative position in x-y space).
It has the shape (number of nodes, maximum number of links at node).
Recall that the standard landlab link direction goes from the tail node to the head node.
A value of zero indicates that the link does not exist or is not active.
A value of -1 indicates that water flow based on
flow__link_to_receiver_node goes from head node to tail node, while
a value of 1 indicates that water flow goes from tail node to head
node.
Examples
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorSteepest
>>> mg = RasterModelGrid((3,3))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fd.run_one_step()
>>> fd.flow_link_direction_at_node()
array([[ 0, 0, 0, 0],
[ 0, -1, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, -1],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0]], dtype=int8)
This method will be updated when the DepressionFinderAndRouter is run.
First, without DepressionFinderAndRouter:
>>> from landlab.components import FlowAccumulator
>>> mg1 = RasterModelGrid((5, 5))
>>> z1 = mg1.add_field(
... "topographic__elevation",
... mg1.x_of_node+2 * mg1.y_of_node,
... at="node",
... )
>>> z1[12] -= 5
>>> mg1.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fa1 = FlowAccumulator(mg1, flow_director='Steepest')
>>> fa1.run_one_step()
>>> fa1.flow_director.links_to_receiver
array([-1, -1, -1, -1, -1,
-1, 5, 15, 7, -1,
-1, 19, -1, 20, -1,
-1, 23, 24, 25, -1,
-1, -1, -1, -1, -1])
>>> fa1.flow_director.flow_link_direction_at_node()
array([[ 0, 0, 0, 0],
[ 0, -1, 0, 0],
[ 0, 0, 0, 0],
[ 0, -1, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, -1],
[ 0, 1, 0, 0],
[ 0, 0, 0, -1],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 1, -1, 0, 0],
[-1, -1, 1, 1],
[ 0, -1, -1, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, -1],
[ 0, 0, 0, -1],
[ 0, 0, 0, -1],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0]], dtype=int8)
Next with DepressionFinderAndRouter:
>>> mg2 = RasterModelGrid((5, 5))
>>> z2 = mg2.add_field(
... "topographic__elevation",
... mg2.x_of_node+2 * mg2.y_of_node,
... at="node",
... )
>>> z2[12] -= 5
>>> mg2.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fa2 = FlowAccumulator(
... mg2,
... flow_director='Steepest',
... depression_finder='DepressionFinderAndRouter',
... routing='D4'
... )
>>> fa2.run_one_step()
>>> fa2.flow_director.links_to_receiver
array([-1, -1, -1, -1, -1,
-1, 5, 6, 7, -1,
-1, 19, 15, 20, -1,
-1, 23, 24, 25, -1,
-1, -1, -1, -1, -1])
>>> fa2.flow_director.flow_link_direction_at_node()
array([[ 0, 0, 0, 0],
[ 0, -1, 0, 0],
[ 0, -1, 0, 0],
[ 0, -1, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, -1],
[ 0, -1, 0, -1],
[ 0, 0, 0, -1],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 1, -1, 0, 0],
[-1, -1, 1, -1],
[ 0, -1, -1, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, -1],
[ 0, 0, 0, -1],
[ 0, 0, 0, -1],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0]], dtype=int8)
flow_link_incoming_at_node()[source]¶Return array that mirrors links at node and indicates incoming flow.
This array has the shape (number of nodes, maximum number of links at node).
Incoming flow is indicated as 1 and outgoing as -1. 0 indicates that no flow moves along the link or that the link does not exist.
Examples
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorSteepest
>>> mg = RasterModelGrid((3,3))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fd.run_one_step()
>>> fd.flow_link_incoming_at_node()
array([[ 0, 0, 0, 0],
[ 0, 1, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, -1],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 0, 0, 0]], dtype=int8)
run_one_step()[source]¶Find flow directions and save to the model grid.
run_one_step() checks for updated boundary conditions, calculates slopes on links, finds baselevel nodes based on the status at node, calculates flow directions, and saves results to the grid.
An alternative to direct_flow() is run_one_step() which does the same things but also returns the receiver nodes not return values.
updated_boundary_conditions()[source]¶Method to update FlowDirectorSteepest when boundary conditions change.
Call this if boundary conditions on the grid are updated after the component is instantiated.
upstream_node_at_link()[source]¶At-link array of the upstream node based on flow direction.
BAD_INDEX_VALUE is given if no upstream node is defined.
Examples
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorSteepest
>>> mg = RasterModelGrid((3,3))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fd.run_one_step()
>>> fd.upstream_node_at_link()
array([-1, -1, -1, 4, -1, -1, -1, -1, -1, -1, -1, -1])
FractureGridGenerator(grid, frac_spacing=10.0, seed=0)[source]¶Bases: landlab.core.model_component.Component
Create a 2D grid with randomly generated fractures.
The grid contains the value 1 where fractures (one cell wide) exist, and 0 elsewhere. The idea is to use this for simulations based on weathering and erosion of, and/or flow within, fracture networks.
Examples
>>> from landlab import RasterModelGrid
>>> grid = RasterModelGrid((5, 5))
>>> fg = FractureGridGenerator(grid=grid, frac_spacing=3)
>>> fg.run_one_step()
>>> grid.at_node['fracture_at_node']
array([1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0], dtype=int8)
Notes
Potential improvements:
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
None Listed
Initialize the FractureGridGenerator.
| Parameters: |
|
|---|
gFlex(grid, Youngs_modulus=650000000000.0, Poissons_ratio=0.25, rho_mantle=3300.0, rho_fill=0.0, elastic_thickness=35000.0, Method='FD', Solver='direct', PlateSolutionType='vWC1994', quiet=True, BC_W='0Displacement0Slope', BC_E='0Displacement0Slope', BC_N='0Displacement0Slope', BC_S='0Displacement0Slope', g=9.80665)[source]¶Bases: landlab.core.model_component.Component
This is a Landlab wrapper for A Wickert’s gFlex flexure model (Wickert et al., 2016, Geoscientific Model Development). The most up-to-date version of his code can be found at github.com/awickert/gFlex.
This Landlab wrapper will use a snapshot of that code, which YOU need to install on your own machine. A stable snapshot of gFlex is hosted on PyPI, which is the recommended version to install. If you have pip (the Python package install tool), simply run ‘pip install gFlex’ from a command prompt. Alternatively, you can download and unpack the code (from github, or with PyPI, pypi.python.org/pypi/gFlex/), then run ‘python setup.py install’.
Note that gFlex maintains its own internal version if the grid, but this should not affect performance.
This component will modify the topographic__elevation field only if one already exists. Note that the gFlex component demands lengths in meters, including the grid dimensions. The component also recognises the gFlex specific parameters ‘Method’, ‘PlateSolutionType’, ‘Solver’, and ‘Quiet’. See the gFlex software documentation for more details.
Examples
NB: these tests are not actually run as our automated testing becomes confused if gFlex is not installed on the testing machine!
>>> from landlab import RasterModelGrid
>>> from landlab.components import gFlex
>>> mg = RasterModelGrid((10, 10), xy_spacing=25000.)
>>> z = mg.add_zeros('topographic__elevation', at='node', dtype=float)
>>> stress = mg.add_zeros('surface_load__stress', at='node', dtype=float)
>>> stress.view().reshape(mg.shape)[3:7, 3:7] += 1.e6
>>> gf = gFlex(mg, BC_E='0Moment0Shear', BC_N='Periodic',
... BC_S='Periodic') # doctest: +SKIP
>>> gf.run_one_step() # doctest: +SKIP
N-S profile across flexed plate:
>>> z.reshape(mg.shape)[:, 5] # doctest: +SKIP
array([-4.54872677, -4.6484927 , -4.82638669, -5.03001546, -5.15351385,
-5.15351385, -5.03001546, -4.82638669, -4.6484927 , -4.54872677])
W-E profile, noting the free BC to the east side:
>>> z.reshape(mg.shape)[5, :] # doctest: +SKIP
array([-0.43536739, -1.19197738, -2.164915 , -3.2388464 , -4.2607558 ,
-5.15351385, -5.89373366, -6.50676947, -7.07880156, -7.63302576])
References
Required Software Citation(s) Specific to this Component
Wickert, A. (2016). Open-source modular solutions for flexural isostasy: gFlex v1.0. Geoscientific Model Development 9(3), 997-1017. https://dx.doi.org/10.5194/gmd-9-997-2016
Additional References
None Listed
Constructor for Wickert’s gFlex in Landlab.
| Parameters: |
|
|---|
__init__(grid, Youngs_modulus=650000000000.0, Poissons_ratio=0.25, rho_mantle=3300.0, rho_fill=0.0, elastic_thickness=35000.0, Method='FD', Solver='direct', PlateSolutionType='vWC1994', quiet=True, BC_W='0Displacement0Slope', BC_E='0Displacement0Slope', BC_N='0Displacement0Slope', BC_S='0Displacement0Slope', g=9.80665)[source]¶Constructor for Wickert’s gFlex in Landlab.
| Parameters: |
|
|---|
GroundwaterDupuitPercolator(grid, hydraulic_conductivity=0.001, porosity=0.2, recharge_rate=1e-08, regularization_f=0.01, courant_coefficient=0.5, vn_coefficient=0.8, callback_fun=<function GroundwaterDupuitPercolator.<lambda>>, **callback_kwds)[source]¶Bases: landlab.core.model_component.Component
Simulate groundwater flow in a shallow unconfined aquifer.
The GroundwaterDupuitPercolator solves the Boussinesq equation for flow in an unconfined aquifer over an impermeable aquifer base and calculates groundwater return flow to the surface. This method uses the Dupuit-Forcheimer approximation. This means that the model assumes the aquifer is laterally extensive in comparison to its thickness, such that the vertical component of flow is negligible. It also assumes that the capillary fringe is small, such that the water table can be modeled as a free surface. Please consider the applicability of these assumptions when using this model. For more details, see component documentation here.
Examples
Import the grid class and component
>>> from landlab import RasterModelGrid
>>> from landlab.components import GroundwaterDupuitPercolator
Initialize the grid and component
>>> grid = RasterModelGrid((10, 10), xy_spacing=10.0)
>>> elev = grid.add_zeros("topographic__elevation", at="node")
>>> abe = grid.add_zeros("aquifer_base__elevation", at="node")
>>> elev[:] = 5.0
>>> gdp = GroundwaterDupuitPercolator(grid)
Run component forward.
>>> dt = 1E4
>>> for i in range(100):
... gdp.run_one_step(dt)
When the model generates groundwater return flow, the surface water flux out of the domain can be calculated only after a FlowAccumulator is run. Below is a more complex model that demonstrates this case.
>>> from landlab.components import FlowAccumulator
Set boundary conditions and initialize grid
>>> grid = RasterModelGrid((5, 41), xy_spacing=10.0)
>>> grid.set_closed_boundaries_at_grid_edges(True, True, False, True)
Make a sloping, 3 m thick aquifer, initially fully saturated
>>> elev = grid.add_zeros("topographic__elevation", at="node")
>>> elev[:] = grid.x_of_node/100+3
>>> base = grid.add_zeros("aquifer_base__elevation", at="node")
>>> base[:] = grid.x_of_node/100
>>> wt = grid.add_zeros("water_table__elevation", at="node")
>>> wt[:] = grid.x_of_node/100 + 3
Initialize components
>>> gdp = GroundwaterDupuitPercolator(grid, recharge_rate=1E-7)
>>> fa = FlowAccumulator(grid, runoff_rate='surface_water__specific_discharge')
Advance timestep. Default units are meters and seconds, though the component is unit agnostic.
>>> dt = 1E3
>>> for i in range(1000):
... gdp.run_one_step(dt)
Calculate surface water flux out of domain
>>> fa.run_one_step()
>>> np.testing.assert_almost_equal(gdp.calc_sw_flux_out(),0.0005077)
Notes
Below is a summary of the theory and numerical implementation of
the GroundwaterDupuitPercolator. A complete description can be found
here.
Groundwater discharge per unit length, \(q\), is calculated as:
where \(K_{sat}\) is the saturated hydraulic conductivity, \(h\) is the aquifer thickness, and \(\alpha\) is the slope angle of the aquifer base.
Surface water discharge per unit area, \(q_s\), is calculated as:
where \(\mathcal{G}_r\) is a smoothed step function, \(\mathcal{R}\) is the ramp function, \(d\) is the regolith thickness, and \(f\) is the recharge rate.
The evolution of aquifer thickness is then given by:
where \(n\) is the drainable porosity.
An explicit forward in time finite volume method is used to implement a numerical solution. Groundwater flow between neighboring nodes is calculated using the saturated thickness at the up-gradient node.
References
Required Software Citation(s) Specific to this Component
Litwin, D. G., Tucker, G.E., Barnhart, K. R., Harman, C. J. (2020). GroundwaterDupuitPercolator: A Landlab component for groundwater flow. Journal of Open Source Software, 5(46), 1935, https://doi.org/10.21105/joss.01935.
Additional References
Marçais, J., de Dreuzy, J. R. & Erhel, J. Dynamic coupling of subsurface and seepage flows solved within a regularized partition formulation. Advances in Water Resources 109, 94–105 (2017).
Childs, E. C. Drainage of Groundwater Resting on a Sloping Bed. Water Resources Research 7, 1256–1263 (1971).
| Parameters: |
|
|---|
K¶hydraulic conductivity at link (m/s)
__init__(grid, hydraulic_conductivity=0.001, porosity=0.2, recharge_rate=1e-08, regularization_f=0.01, courant_coefficient=0.5, vn_coefficient=0.8, callback_fun=<function GroundwaterDupuitPercolator.<lambda>>, **callback_kwds)[source]¶| Parameters: |
|
|---|
calc_gw_flux_out()[source]¶Groundwater flux through open boundaries may be positive (out of the domain) or negative (into the domain).
This function determines the correct sign for specific discharge based upon this convention, and sums the flux across the boundary faces. (m3/s)
calc_recharge_flux_in()[source]¶Calculate flux into the domain from recharge.
Includes recharge that may immediately become saturation excess overland flow. (m3/s)
calc_sw_flux_out()[source]¶Surface water flux out of the domain through seepage and saturation excess.
Note that model does not allow for reinfiltration. (m3/s)
callback_fun¶callback function for adaptive timestep solver
| Parameters: | callback_fun (function(grid, recharge_rate, substep_dt, **callback_kwds)) – Optional function that will be executed at the end of each sub-timestep in the run_with_adaptive_time_step_solver method. Intended purpose is to write output not otherwise visible outside of the method call. The function should have three required arguments: grid: the ModelGrid instance used by GroundwaterDupuitPercolator recharge_rate: an array at node that is the specified recharge rate substep_dt: the length of the current substep determined internally by run_with_adaptive_time_step_solver to meet stability criteria. Callback functions with two arguments (grid, substep_dt) are depricated. |
|---|
courant_coefficient¶Courant coefficient for adaptive time step.
| Parameters: | courant_coefficient (float (-)) – The muliplying factor on the condition that the timestep is
smaller than the minimum link length over groundwater flow
velocity. This parameter is only used with
run_with_adaptive_time_step_solver and must be greater than
zero. |
|---|
n¶drainable porosity of the aquifer (-)
number_of_substeps¶The number of substeps used by the run_with_adaptive_time_step_solver method in the latest method call.
recharge¶recharge rate (m/s)
run_one_step(dt)[source]¶Advance component by one time step of size dt.
| Parameters: | dt (float) – The imposed timestep. |
|---|
run_with_adaptive_time_step_solver(dt)[source]¶Advance component by one time step of size dt, subdividing the timestep into substeps as necessary to meet stability conditions. Note this method returns the fluxes at the last substep, but also returns a new field, average_surface_water__specific_discharge, that is averaged over all subtimesteps. To return state during substeps, provide a callback_fun.
| Parameters: | dt (float) – The imposed timestep. |
|---|
vn_coefficient¶Coefficient for the diffusive timestep condition in the adaptive timestep solver.
| Parameters: | vn_coefficient (float (-)) – The multiplying factor C for the condition dt >= C*dx^2/(4D),
where D = Kh/n is the diffusivity of the Boussinesq
equation. This arises from a von Neumann stability analysis of
the Boussinesq equation when the hydraulic gradient is small.
This parameter is only used with run_with_adaptive_time_step_solver
and must be greater than zero. |
|---|
HackCalculator(grid, save_full_df=False, **kwds)[source]¶Bases: landlab.core.model_component.Component
This component calculates Hack’s law parameters for drainage basins.
Hacks law is given as
Where \(L\) is the distance to the drainage divide along the channel, \(A\) is the drainage area, and \(C\) are parameters.
The HackCalculator uses a ChannelProfiler to determine the nodes on which to calculate the parameter fit.
Examples
>>> import pandas as pd
>>> pd.set_option('display.max_columns', None)
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import (
... FlowAccumulator,
... FastscapeEroder,
... HackCalculator)
>>> np.random.seed(42)
>>> mg = RasterModelGrid((50, 100), xy_spacing=100)
>>> z = mg.add_zeros('node', 'topographic__elevation')
>>> z[mg.core_nodes] += np.random.randn(mg.core_nodes.size)
>>> fa = FlowAccumulator(mg)
>>> fs = FastscapeEroder(mg, K_sp=0.001)
>>> for i in range(100):
... fa.run_one_step()
... fs.run_one_step(1000)
... z[mg.core_nodes] += 0.01 * 1000
>>> hc = HackCalculator(mg)
>>> hc.calculate_hack_parameters()
>>> largest_outlet = mg.boundary_nodes[
... np.argsort(mg.at_node['drainage_area'][mg.boundary_nodes])[-1:]][0]
>>> largest_outlet
4978
>>> hc.hack_coefficient_dataframe.loc[largest_outlet, "A_max"]
2830000.0
>>> hc.hack_coefficient_dataframe.round(2) # doctest: +NORMALIZE_WHITESPACE
A_max C h
basin_outlet_id
4978 2830000.0 0.31 0.62
>>> hc = HackCalculator(
... mg,
... number_of_watersheds=3,
... main_channel_only=False,
... save_full_df=True)
>>> hc.calculate_hack_parameters()
>>> hc.hack_coefficient_dataframe.round(2) # doctest: +NORMALIZE_WHITESPACE
A_max C h
basin_outlet_id
39 2170000.0 0.13 0.69
4929 2350000.0 0.13 0.68
4978 2830000.0 0.23 0.64
>>> hc.full_hack_dataframe.head().round(2) # doctest: +NORMALIZE_WHITESPACE
basin_outlet_id A L_obs L_est
node_id
39 39.0 2170000.0 3200.0 2903.43
139 39.0 2170000.0 3100.0 2903.43
238 39.0 10000.0 0.0 71.61
239 39.0 2160000.0 3000.0 2894.22
240 39.0 10000.0 0.0 71.61
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
Hack, J. T. Studies of longitudinal stream profiles in Virginia and Maryland (Vol. 294). U.S. Geological Survey Professional Paper 294-B (1957). https://doi.org/10.3133/pp294B
| Parameters: |
|
|---|
__init__(grid, save_full_df=False, **kwds)[source]¶| Parameters: |
|
|---|
full_hack_dataframe¶Full Hack calculation dataframe.
This dataframe is optionally created and stored on the component when
the keyword argument full_hack_dataframe=True is passed to the
component init.
It is pandas dataframe with a row for every model grid cell used to estimate the Hack parameters. It has the following index and columns.
- Index
- node_id*: The node ID of the model grid cell.
- Columns
- basin_outlet_id: The node IDs of watershed outlet
- A: The drainage are of the model grid cell.
- L_obs: The observed distance to the divide.
- L_est: The predicted distance to divide based on the Hack coefficient fit.
hack_coefficient_dataframe¶Hack coefficient dataframe.
This dataframe is created and stored on the component.
It is a pandas dataframe with one row for each basin for which Hack parameters are calculated. Thus, there are as many rows as the number of watersheds identified by the ChannelProfiler.
The dataframe has the following index and columns.
- Index
- basin_outlet_id: The node ID of the watershed outlet where each set of Hack parameters was estimated.
- Columns
- A_max: The drainage area of the watershed outlet.
- C: The Hack coefficient as defined in the equations above.
- h: The Hack exponent as defined in the equations above.
HeightAboveDrainageCalculator(grid, channel_mask='channel__mask')[source]¶Bases: landlab.core.model_component.Component
Calculate the elevation difference between each node and its nearest drainage node in a DEM.
This component implements the method described by Nobre et al (2011). A single direction flow director (D8 or steepest descent) must be run prior to HeightAboveDrainageCalculator to supply the flow directions. This component does not fill depressions in a DEM, but rather it treats them as drainage nodes. For best results, please run one of the available pit filling components prior to HeightAboveDrainageCalculator.
Examples
>>> import numpy as np
>>> from numpy.testing import assert_equal
>>> from landlab import RasterModelGrid
>>> from landlab.components import HeightAboveDrainageCalculator, FlowAccumulator
>>> mg = RasterModelGrid((4, 5))
>>> z = mg.add_zeros("topographic__elevation", at="node")
>>> mg.set_status_at_node_on_edges(right=mg.BC_NODE_IS_CLOSED, bottom=mg.BC_NODE_IS_FIXED_VALUE, left=mg.BC_NODE_IS_CLOSED, top=mg.BC_NODE_IS_CLOSED)
>>> elev = np.array([[2,1,0,1,2],[3,2,1,2,3],[4,3,2,3,4],[5,4,4,4,5]])
>>> z[:] = elev.reshape(len(z))
>>> elev
array([[2, 1, 0, 1, 2],
[3, 2, 1, 2, 3],
[4, 3, 2, 3, 4],
[5, 4, 4, 4, 5]])
>>> fa = FlowAccumulator(mg, flow_director="D8")
>>> fa.run_one_step()
>>> channel__mask = mg.zeros(at="node")
>>> channel__mask[[2,7]] = 1
>>> channel__mask.reshape(elev.shape)
array([[ 0., 0., 1., 0., 0.],
[ 0., 0., 1., 0., 0.],
[ 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0.]])
>>> hd = HeightAboveDrainageCalculator(mg, channel_mask=channel__mask)
>>> hd.run_one_step()
>>> mg.at_node["height_above_drainage__elevation"].reshape(elev.shape) # doctest: +NORMALIZE_WHITESPACE
array([[ 2., 0., 0., 0., 0.],
[ 3., 2., 0., 2., 3.],
[ 4., 2., 1., 2., 4.],
[ 5., 4., 4., 4., 5.]])
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
Nobre, A. D., Cuartas, L. A., Hodnett, M., Rennó, C. D., Rodrigues, G., Silveira, A., et al. (2011). Height Above the Nearest Drainage – a hydrologically relevant new terrain model. Journal of Hydrology, 404(1), 13–29. https://doi.org/10.1016/j.jhydrol.2011.03.051
| Parameters: |
|
|---|
__init__(grid, channel_mask='channel__mask')[source]¶| Parameters: |
|
|---|
channel_mask¶KinwaveImplicitOverlandFlow(grid, runoff_rate=1.0, roughness=0.01, changing_topo=False, depth_exp=1.5, weight=1.0)[source]¶Bases: landlab.core.model_component.Component
Calculate shallow water flow over topography.
Landlab component that implements a two-dimensional kinematic wave model. This is a form of the 2D shallow-water equations in which energy slope is assumed to equal bed slope. The solution method is locally implicit, and works as follows. At each time step, we iterate from upstream to downstream over the topography. Because we are working downstream, we can assume that we know the total water inflow to a given cell. We solve the following mass conservation equation at each cell:
where \(H\) is water depth, \(t\) indicates time step number, \(\Delta t\) is time step duration, \(Q_{in}\) is total inflow discharge, \(Q_{out}\) is total outflow discharge, \(A\) is cell area, and \(R\) is local runoff rate (precipitation minus infiltration; could be negative if runon infiltration is occurring).
The specific outflow discharge leaving a cell along one of its faces is:
where \(C_r\) is a roughness coefficient (such as Manning’s n), \(\alpha\) is an exponent equal to \(5/3\) for the Manning equation and \(3/2\) for the Chezy family, and \(S\) is the downhill-positive gradient of the link that crosses this particular face. Outflow discharge is zero for links that are flat or “uphill” from the given node. Total discharge out of a cell is then the sum of (specific discharge x face width) over all outflow faces
where \(N\) is the number of outflow faces (i.e., faces where the ground slopes downhill away from the cell’s node), and \(W_i\) is the width of face \(i\).
We use the depth at the cell’s node, so this simplifies to:
We define \(H\) in the above as a weighted sum of the “old” (time step \(t\)) and “new” (time step \(t+1\)) depth values:
If \(w=1\), the method is fully implicit. If \(w=0\), it is a simple forward explicit method.
When we combine these equations, we have an equation that includes the unknown \(H^{t+1}\) and a bunch of terms that are known. If \(w\ne 0\), it is a nonlinear equation in \(H^{t+1}\), and must be solved iteratively. We do this using a root-finding method in the scipy.optimize library.
Examples
>>> from landlab import RasterModelGrid
>>> rg = RasterModelGrid((4, 5), xy_spacing=10.0)
>>> z = rg.add_zeros("topographic__elevation", at="node")
>>> kw = KinwaveImplicitOverlandFlow(rg)
>>> round(kw.runoff_rate * 1.0e7, 2)
2.78
>>> kw.vel_coef # default value
100.0
>>> rg.at_node['surface_water__depth'][6:9]
array([ 0., 0., 0.])
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
None Listed
Initialize the KinwaveImplicitOverlandFlow.
| Parameters: |
|
|---|
__init__(grid, runoff_rate=1.0, roughness=0.01, changing_topo=False, depth_exp=1.5, weight=1.0)[source]¶Initialize the KinwaveImplicitOverlandFlow.
| Parameters: |
|
|---|
depth¶The depth of water at each node.
runoff_rate¶Runoff rate.
| Parameters: | runoff_rate (float, optional (defaults to 1 mm/hr)) – Precipitation rate, mm/hr. The value provide is divided by 3600000.0. |
|---|---|
| Returns: | |
| Return type: | The current value of the runoff rate. |
vel_coef¶Velocity coefficient.
KinwaveOverlandFlowModel(grid, precip_rate=1.0, precip_duration=1.0, infilt_rate=0.0, roughness=0.01)[source]¶Bases: landlab.core.model_component.Component
Calculate water flow over topography.
Landlab component that implements a two-dimensional kinematic wave model. This is an extremely simple, unsophisticated model, originally built simply to demonstrate the component creation process. Limitations to the present version include: infiltration is handled very crudely, the called is responsible for picking a stable time step size (no adaptive time stepping is used in the run_one_step method), precipitation rate is constant for a given duration (then zero), and all parameters are uniform in space. Also, the terrain is assumed to be stable over time. Caveat emptor!
Examples
>>> from landlab import RasterModelGrid
>>> rg = RasterModelGrid((4, 5), xy_spacing=10.0)
>>> z = rg.add_zeros("topographic__elevation", at="node")
>>> s = rg.add_zeros("topographic__gradient", at="link")
>>> kw = KinwaveOverlandFlowModel(rg)
>>> kw.vel_coef
100.0
>>> rg.at_node['surface_water__depth']
array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0.])
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
None Listed
Initialize the KinwaveOverlandFlowModel.
| Parameters: |
|
|---|
__init__(grid, precip_rate=1.0, precip_duration=1.0, infilt_rate=0.0, roughness=0.01)[source]¶Initialize the KinwaveOverlandFlowModel.
| Parameters: |
|
|---|
run_one_step(dt)[source]¶Calculate water flow for a time period dt.
Default units for dt are seconds.
vel_coef¶Velocity coefficient.
(1/roughness)
LakeMapperBarnes(grid, surface='topographic__elevation', method='Steepest', fill_flat=True, fill_surface='topographic__elevation', redirect_flow_steepest_descent=False, reaccumulate_flow=False, ignore_overfill=False, track_lakes=True)[source]¶Bases: landlab.core.model_component.Component
A Landlab implementation of the Barnes et al. (2014) lake filling & lake routing algorithms, lightly modified and adapted for Landlab by DEJH. This component is designed as a direct replacement for the LakeMapper as existing pre-Aug 2018, and provides a suite of properties to access information about the lakes created each time it is run. Only significant difference is the way the lakes are coded: this component uses the (unique) ID of the outlet node, whereas DepressionFinderAndRouter uses one of the pit node IDs. Note also this component does not offer the lake_codes or display_depression_map options, for essentially this reason. Use lake_map instead for both. It also uses a much more Landlabbian run_one_step() method as its driver, superceding DepressionFinderAndRouter’s map_depressions().
A variety of options is provided. Flow routing is route-to-one in this implementation, but can be either D4 (“steepest”) or D8 on a raster. The surface can be filled to either flat or a very slight downward incline, such that subsequent flow routing will run over the lake surface. This incline is applied at machine precision to minimise the chances of creating false outlets by overfill, and note that the gradient as calculated on such surfaces may still appear to be zero. The filling can either be performed in place, or on a new (water) surface distinct from the original (rock) surface. For efficiency, data structures describing the lakes and their properties are only created, and existing flow direction and flow accumulation fields modified, if those flags are set at instantiation.
With care, this component can be used to create a dynamically flooding surface in a fluvial landscape (interacting with, e.g., the StreamPowerEroder). See the run_one_step docstring for an example.
References
Required Software Citation(s) Specific to this Component
Barnes, R., Lehman, C., Mulla, D. (2014). Priority-flood: An optimal depression-filling and watershed-labeling algorithm for digital elevation models. Computers and Geosciences 62(C), 117 - 127. https://dx.doi.org/10.1016/j.cageo.2013.04.024
Additional References
None Listed
Initialize the component.
| Parameters: |
|
|---|
__init__(grid, surface='topographic__elevation', method='Steepest', fill_flat=True, fill_surface='topographic__elevation', redirect_flow_steepest_descent=False, reaccumulate_flow=False, ignore_overfill=False, track_lakes=True)[source]¶Initialize the component.
| Parameters: |
|
|---|
lake_areas¶A nlakes-long array of the area of each lake. The order is the same as that of the keys in lake_dict, and of lake_outlets. Note that outlet nodes are not parts of the lakes.
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import LakeMapperBarnes, FlowAccumulator
>>> mg = RasterModelGrid((5, 6))
>>> for edge in ('left', 'top', 'bottom'):
... mg.status_at_node[mg.nodes_at_edge(edge)] = mg.BC_NODE_IS_CLOSED
>>> z = mg.add_zeros("topographic__elevation", at="node", dtype=float)
>>> z[:] = mg.node_x.max() - mg.node_x
>>> z[[10, 23]] = 1.1 # raise "guard" exit nodes
>>> z[7] = 2. # is a lake on its own
>>> z[9] = 0.5
>>> z[15] = 0.3
>>> z[14] = 0.6 # [9, 14, 15] is a lake
>>> z[22] = 0.9 # a non-contiguous lake node also draining to 16
>>> z_init = z.copy()
>>> fa = FlowAccumulator(mg)
>>> lmb = LakeMapperBarnes(mg, method='Steepest', fill_flat=False,
... redirect_flow_steepest_descent=False,
... track_lakes=False)
>>> lmb.run_one_step()
>>> try:
... lmb.lake_areas # note track_lakes=False
... except ValueError:
... print('ValueError was raised: ' +
... 'Enable tracking to access information about lakes')
ValueError was raised: Enable tracking to access information about lakes
>>> z[:] = z_init
>>> lmb = LakeMapperBarnes(mg, method='Steepest', fill_flat=False,
... redirect_flow_steepest_descent=False,
... track_lakes=True)
>>> lmb.run_one_step()
>>> lmb.lake_outlets
[16, 8]
>>> np.allclose(lmb.lake_areas, np.array([ 4., 1.]))
True
lake_at_node¶Return a boolean array, True if the node is flooded, False otherwise. The outlet nodes are NOT part of the lakes.
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import LakeMapperBarnes, FlowAccumulator
>>> mg = RasterModelGrid((5, 6))
>>> for edge in ('left', 'top', 'bottom'):
... mg.status_at_node[mg.nodes_at_edge(edge)] = mg.BC_NODE_IS_CLOSED
>>> z = mg.add_zeros("topographic__elevation", at="node", dtype=float)
>>> z[:] = mg.node_x.max() - mg.node_x
>>> z[[10, 23]] = 1.1 # raise "guard" exit nodes
>>> z[7] = 2. # is a lake on its own
>>> z[9] = 0.5
>>> z[15] = 0.3
>>> z[14] = 0.6 # [9, 14, 15] is a lake
>>> z[22] = 0.9 # a non-contiguous lake node also draining to 16
>>> z_init = z.copy()
>>> fa = FlowAccumulator(mg)
>>> lmb = LakeMapperBarnes(mg, method='Steepest', fill_flat=False,
... redirect_flow_steepest_descent=False,
... track_lakes=True)
>>> lmb.run_one_step()
>>> lake_at_node = np.array([False, False, False, False, False, False,
... False, True, False, True, False, False,
... False, False, True, True, False, False,
... False, False, False, False, True, False,
... False, False, False, False, False, False],
... dtype=bool)
>>> np.all(np.equal(lmb.lake_at_node, lake_at_node))
True
lake_depths¶Return the change in surface elevation at each node this step. Requires that surface and fill_surface were not the same array at instantiation.
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import LakeMapperBarnes, FlowAccumulator
>>> mg = RasterModelGrid((5, 6))
>>> for edge in ('left', 'top', 'bottom'):
... mg.status_at_node[mg.nodes_at_edge(edge)] = mg.BC_NODE_IS_CLOSED
>>> z = mg.add_zeros("topographic__elevation", at="node", dtype=float)
>>> z[:] = mg.node_x.max() - mg.node_x
>>> z[[10, 23]] = 1.1 # raise "guard" exit nodes
>>> z[7] = 2. # is a lake on its own
>>> z[9] = 0.5
>>> z[15] = 0.3
>>> z[14] = 0.6 # [9, 14, 15] is a lake
>>> z[22] = 0.9 # a non-contiguous lake node also draining to 16
>>> z_init = z.copy()
>>> fa = FlowAccumulator(mg)
>>> lmb = LakeMapperBarnes(mg, method='Steepest', fill_flat=False,
... redirect_flow_steepest_descent=False,
... track_lakes=True)
>>> lmb.run_one_step()
>>> try: # won't work as surface & fill_surface are both z
... lmb.lake_depths
... except ValueError:
... print('ValueError was raised: ' +
... 'surface and fill_surface must be different fields ' +
... 'or arrays to enable the property lake_depths!')
ValueError was raised: surface and fill_surface must be different fields or arrays to enable the property lake_depths!
>>> z[:] = z_init
>>> lmb = LakeMapperBarnes(mg, method='Steepest', fill_flat=False,
... surface=z_init,
... redirect_flow_steepest_descent=False,
... track_lakes=True)
>>> lmb.run_one_step()
>>> lake_depths = np.array([ 0. , 0. , 0. , 0. , 0. , 0. ,
... 0. , 1. , 0. , 0.5, 0. , 0. ,
... 0. , 0. , 0.4, 0.7, 0. , 0. ,
... 0. , 0. , 0. , 0. , 0.1, 0. ,
... 0. , 0. , 0. , 0. , 0. , 0. ])
>>> np.all(np.equal(lmb.lake_depths,
... lake_depths)) # slope applied, so...
False
>>> np.allclose(lmb.lake_depths, lake_depths)
True
lake_dict¶Return a dictionary where the keys are the outlet nodes of each lake, and the values are deques of nodes within each lake. Items are not returned in ID order. The outlet nodes are NOT part of the lake.
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import LakeMapperBarnes, FlowAccumulator
>>> mg = RasterModelGrid((5, 6))
>>> for edge in ('left', 'top', 'bottom'):
... mg.status_at_node[mg.nodes_at_edge(edge)] = mg.BC_NODE_IS_CLOSED
>>> z = mg.add_zeros("topographic__elevation", at="node", dtype=float)
>>> z.reshape(mg.shape)[2, 1:-1] = [2., 1., 0.5, 1.5]
>>> z.reshape(mg.shape)[1, 1:-1] = [2.1, 1.1, 0.6, 1.6]
>>> z.reshape(mg.shape)[3, 1:-1] = [2.2, 1.2, 0.7, 1.7]
>>> z_init = z.copy()
>>> fa = FlowAccumulator(mg)
>>> lmb = LakeMapperBarnes(mg, method='Steepest', surface=z_init,
... fill_surface=z, fill_flat=False,
... redirect_flow_steepest_descent=False,
... track_lakes=False)
>>> lmb.run_one_step()
>>> try:
... lmb.lake_dict
... except ValueError:
... print('ValueError was raised: ' +
... 'Enable tracking to access information about lakes')
ValueError was raised: Enable tracking to access information about lakes
>>> lmb = LakeMapperBarnes(mg, method='Steepest', surface=z_init,
... fill_surface=z, fill_flat=False,
... redirect_flow_steepest_descent=False,
... track_lakes=True)
>>> lmb.run_one_step()
>>> lmb.lake_dict == {16: deque([15, 9, 8, 14, 20, 21])}
True
lake_map¶Return an array of ints, where each node within a lake is labelled with its outlet node ID. The outlet nodes are NOT part of the lakes. Nodes not in a lake are labelled with BAD_INDEX_VALUE (default -1).
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import LakeMapperBarnes, FlowAccumulator
>>> mg = RasterModelGrid((5, 6))
>>> for edge in ('left', 'top', 'bottom'):
... mg.status_at_node[mg.nodes_at_edge(edge)] = mg.BC_NODE_IS_CLOSED
>>> z = mg.add_zeros("topographic__elevation", at="node", dtype=float)
>>> z[:] = mg.node_x.max() - mg.node_x
>>> z[[10, 23]] = 1.1 # raise "guard" exit nodes
>>> z[7] = 2. # is a lake on its own
>>> z[9] = 0.5
>>> z[15] = 0.3
>>> z[14] = 0.6 # [9, 14, 15] is a lake
>>> z[22] = 0.9 # a non-contiguous lake node also draining to 16
>>> z_init = z.copy()
>>> fa = FlowAccumulator(mg)
>>> lmb = LakeMapperBarnes(mg, method='Steepest', fill_flat=False,
... redirect_flow_steepest_descent=False,
... track_lakes=False)
>>> lmb.run_one_step()
>>> try:
... lmb.lake_map
... except ValueError:
... print('ValueError was raised: ' +
... 'Enable tracking to access information about lakes')
ValueError was raised: Enable tracking to access information about lakes
>>> z[:] = z_init
>>> lmb = LakeMapperBarnes(mg, method='Steepest', fill_flat=False,
... redirect_flow_steepest_descent=False,
... track_lakes=True)
>>> lmb.run_one_step()
>>> lake_map = np.array([-1, -1, -1, -1, -1, -1,
... -1, 8, -1, 16, -1, -1,
... -1, -1, 16, 16, -1, -1,
... -1, -1, -1, -1, 16, -1,
... -1, -1, -1, -1, -1, -1])
>>> np.all(np.equal(lmb.lake_map, lake_map))
True
Note that the outlet node is NOT part of the lake.
Updating the elevations works fine:
>>> z.reshape(mg.shape)[2, 1:-1] = [2., 1., 0.5, 1.5]
>>> z.reshape(mg.shape)[1, 1:-1] = [2.1, 1.1, 0.6, 1.6]
>>> z.reshape(mg.shape)[3, 1:-1] = [2.2, 1.2, 0.7, 1.7]
>>> lmb.run_one_step()
>>> new_lake_map = np.array([-1, -1, -1, -1, -1, -1,
... -1, -1, 16, 16, -1, -1,
... -1, -1, 16, 16, -1, -1,
... -1, -1, 16, 16, -1, -1,
... -1, -1, -1, -1, -1, -1])
>>> np.all(np.equal(lmb.lake_map, new_lake_map))
True
lake_outlets¶Returns the outlet for each lake, not necessarily in ID order.
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import LakeMapperBarnes, FlowAccumulator
>>> mg = RasterModelGrid((5, 6))
>>> for edge in ('left', 'top', 'bottom'):
... mg.status_at_node[mg.nodes_at_edge(edge)] = mg.BC_NODE_IS_CLOSED
>>> z = mg.add_zeros("topographic__elevation", at="node", dtype=float)
>>> z.reshape(mg.shape)[2, 1:-1] = [2., 1., 0.5, 1.5]
>>> z.reshape(mg.shape)[1, 1:-1] = [2.1, 1.1, 0.6, 1.6]
>>> z.reshape(mg.shape)[3, 1:-1] = [2.2, 1.2, 0.7, 1.7]
>>> z_init = z.copy()
>>> fa = FlowAccumulator(mg)
>>> lmb = LakeMapperBarnes(mg, method='Steepest', surface=z_init,
... fill_surface=z, fill_flat=False,
... redirect_flow_steepest_descent=False,
... track_lakes=False)
>>> lmb.run_one_step()
>>> try:
... lmb.lake_outlets
... except ValueError:
... print('ValueError was raised: ' +
... 'Enable tracking to access information about lakes')
ValueError was raised: Enable tracking to access information about lakes
>>> lmb = LakeMapperBarnes(mg, method='Steepest', surface=z_init,
... fill_surface=z, fill_flat=False,
... redirect_flow_steepest_descent=False,
... track_lakes=True)
>>> lmb.run_one_step()
>>> lmb.lake_outlets == [16, ]
True
lake_volumes¶A nlakes-long array of the volume of each lake. The order is the same as that of the keys in lake_dict, and of lake_outlets. Note that this calculation is performed relative to the initial surface, so is only a true lake volume if the initial surface was the rock suface (not an earlier water level).
Requires that surface and fill_surface were not the same array at instantiation.
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import LakeMapperBarnes, FlowAccumulator
>>> mg = RasterModelGrid((5, 6))
>>> for edge in ('left', 'top', 'bottom'):
... mg.status_at_node[mg.nodes_at_edge(edge)] = mg.BC_NODE_IS_CLOSED
>>> z = mg.add_zeros("topographic__elevation", at="node", dtype=float)
>>> z[:] = mg.node_x.max() - mg.node_x
>>> z[[10, 23]] = 1.1 # raise "guard" exit nodes
>>> z[7] = 2. # is a lake on its own
>>> z[9] = 0.5
>>> z[15] = 0.3
>>> z[14] = 0.6 # [9, 14, 15] is a lake
>>> z[22] = 0.9 # a non-contiguous lake node also draining to 16
>>> z_init = z.copy()
>>> fa = FlowAccumulator(mg)
>>> lmb = LakeMapperBarnes(mg, method='Steepest', fill_flat=False,
... redirect_flow_steepest_descent=False,
... track_lakes=True)
>>> lmb.run_one_step()
>>> try: # won't work as surface & fill_surface are both z
... lmb.lake_volumes
... except ValueError:
... print('ValueError was raised: ' +
... 'surface and fill_surface must be different fields ' +
... 'or arrays to enable the property lake_volumes!')
ValueError was raised: surface and fill_surface must be different fields or arrays to enable the property lake_volumes!
>>> z[:] = z_init
>>> lmb = LakeMapperBarnes(mg, method='Steepest', fill_flat=False,
... surface=z_init,
... redirect_flow_steepest_descent=False,
... track_lakes=True)
>>> lmb.run_one_step()
>>> lmb.lake_outlets
[16, 8]
>>> np.allclose(lmb.lake_volumes, np.array([ 1.7, 1. ]))
True
number_of_lakes¶Return the number of individual lakes. Lakes sharing outlet nodes are considered part of the same lake.
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import LakeMapperBarnes, FlowAccumulator
>>> mg = RasterModelGrid((5, 6))
>>> for edge in ('left', 'top', 'bottom'):
... mg.status_at_node[mg.nodes_at_edge(edge)] = mg.BC_NODE_IS_CLOSED
>>> z = mg.add_zeros("topographic__elevation", at="node", dtype=float)
>>> z[:] = mg.node_x.max() - mg.node_x
>>> z[[10, 23]] = 1.1 # raise "guard" exit nodes
>>> z[7] = 2. # is a lake on its own
>>> z[9] = 0.5
>>> z[15] = 0.3
>>> z[14] = 0.6 # [9, 14, 15] is a lake
>>> z[22] = 0.9 # a non-contiguous lake node also draining to 16
>>> z_init = z.copy()
>>> fa = FlowAccumulator(mg)
>>> lmb = LakeMapperBarnes(mg, method='Steepest', surface=z_init,
... fill_surface=z, fill_flat=False,
... redirect_flow_steepest_descent=False,
... track_lakes=False)
>>> lmb.run_one_step()
>>> try:
... lmb.number_of_lakes
... except ValueError:
... print('ValueError was raised: ' +
... 'Enable tracking to access information about lakes')
ValueError was raised: Enable tracking to access information about lakes
>>> lmb = LakeMapperBarnes(mg, method='Steepest', surface=z_init,
... fill_surface=z, fill_flat=False,
... redirect_flow_steepest_descent=False,
... track_lakes=True)
>>> lmb.run_one_step()
>>> lmb.number_of_lakes
2
run_one_step()[source]¶Fills the surface to fill all pits. Note that a second run on a surface that has already been filled will not “see” any existing.
lakes correctly - it will see lakes, but with zero depths. In particular, if fill_flat is False, an attempt to fill a surface a second time will raise a ValueError unless ignore_overfill. (In this case, setting ignore_overfill is True will give the expected behaviour.)
If reaccumulate_flow was True at instantiation, run_one_step also updates all the various flow fields produced by the FlowDirector and FlowAccumulator components.
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import LakeMapperBarnes, FlowAccumulator
>>> from landlab.components import FlowDirectorSteepest
>>> mg = RasterModelGrid((5, 6), xy_spacing=2.)
>>> for edge in ('left', 'top', 'bottom'):
... mg.status_at_node[mg.nodes_at_edge(edge)] = mg.BC_NODE_IS_CLOSED
>>> z = mg.add_zeros("topographic__elevation", at="node", dtype=float)
>>> z.reshape(mg.shape)[2, 1:-1] = [2., 1., 0.5, 1.5]
>>> z.reshape(mg.shape)[1, 1:-1] = [2.1, 1.1, 0.6, 1.6]
>>> z.reshape(mg.shape)[3, 1:-1] = [2.2, 1.2, 0.7, 1.7]
>>> z_init = z.copy()
>>> fa = FlowAccumulator(mg)
>>> lmb = LakeMapperBarnes(mg, method='Steepest', surface=z_init,
... fill_flat=False,
... redirect_flow_steepest_descent=False,
... track_lakes=False)
TODO: once return_array_at_node is fixed, this example should also take fill_surface…
>>> lmb.run_one_step()
>>> z_out = np.array([ 0. , 0. , 0. , 0. , 0. , 0. ,
... 0. , 2.1, 1.5, 1.5, 1.6, 0. ,
... 0. , 2. , 1.5, 1.5, 1.5, 0. ,
... 0. , 2.2, 1.5, 1.5, 1.7, 0. ,
... 0. , 0. , 0. , 0. , 0. , 0. ])
>>> np.allclose(z, z_out)
True
>>> np.all(np.equal(z, z_out)) # those 1.5's are actually a bit > 1.5
False
>>> try:
... lmb.lake_map # not created, as we aren't tracking
... except ValueError:
... print('ValueError was raised: ' +
... 'Enable tracking to access information about lakes')
ValueError was raised: Enable tracking to access information about lakes
>>> lmb.was_there_overfill # everything fine with slope adding
False
>>> fd = FlowDirectorSteepest(mg)
>>> fa = FlowAccumulator(mg) # routing will work fine now
>>> fd.run_one_step()
>>> fa.run_one_step()
>>> np.all(mg.at_node['flow__sink_flag'][mg.core_nodes] == 0)
True
>>> drainage_area = np.array([ 0., 0., 0., 0., 0., 0.,
... 0., 4., 8., 12., 4., 4.,
... 0., 4., 16., 36., 40., 40.,
... 0., 4., 8., 4., 4., 4.,
... 0., 0., 0., 0., 0., 0.])
>>> np.allclose(mg.at_node['drainage_area'], drainage_area)
True
Test two pits:
>>> z[:] = mg.node_x.max() - mg.node_x
>>> z[23] = 1.3
>>> z[15] = 0.3
>>> z[10] = 1.3 # raise "guard" exit nodes
>>> z[7] = 2. # is a lake on its own, if D8
>>> z[9] = 0.5
>>> z[14] = 0.6 # [9, 14, 15] is a lake in both methods
>>> z[16] = 1.2
>>> z[22] = 0.9 # a non-contiguous lake node also draining to 16 if D8
>>> z_init = z.copy()
>>> fa = FlowAccumulator(mg)
>>> lmb = LakeMapperBarnes(mg, method='D8', fill_flat=True,
... track_lakes=True)
>>> lmb.run_one_step() # note the D8 routing now
>>> lmb.lake_dict == {22: deque([15, 9, 14])}
True
>>> lmb.number_of_lakes
1
>>> try:
... lmb.lake_depths # z was both surface and 'fill_surface'
... except ValueError:
... print('ValueError was raised: ' +
... 'surface and fill_surface must be different fields ' +
... 'or arrays to enable the property fill_depth!')
ValueError was raised: surface and fill_surface must be different fields or arrays to enable the property fill_depth!
>>> z[:] = z_init
>>> lmb = LakeMapperBarnes(mg, method='Steepest',
... fill_flat=False, track_lakes=True)
>>> lmb.run_one_step() # compare to the method='D8' lakes, above...
>>> lmb.lake_dict == {8: deque([7]), 16: deque([15, 9, 14, 22])}
True
>>> lmb.number_of_lakes
2
>>> np.allclose(lmb.lake_areas, np.array([ 16., 4.]))
True
>>> try:
... lmb.run_one_step() # z already filled, so...
... except ValueError:
... print('ValueError was raised: ' +
... 'Pit is overfilled due to creation of two outlets as ' +
... 'the minimum gradient gets applied. Suppress this ' +
... 'Error with the ignore_overfill flag at component ' +
... 'instantiation.')
ValueError was raised: Pit is overfilled due to creation of two outlets as the minimum gradient gets applied. Suppress this Error with the ignore_overfill flag at component instantiation.
Suppress this behaviour with ignore_overfill:
>>> z[:] = z_init
>>> lmb = LakeMapperBarnes(mg, method='Steepest',
... fill_flat=False, track_lakes=True,
... ignore_overfill=True)
>>> lmb.run_one_step()
>>> lmb.lake_dict == {8: deque([7]), 16: deque([15, 9, 14, 22])}
True
>>> lmb.run_one_step()
>>> np.allclose(lmb.lake_areas, np.array([ 16., 4.])) # found them!
True
The component can redirect flow to account for the fills that have been carried out (all necessary fields get updated):
>>> z[:] = z_init
>>> fd.run_one_step()
>>> init_flowdirs = mg.at_node['flow__receiver_node'].copy()
>>> fa.run_one_step()
>>> init_areas = mg.at_node['drainage_area'].copy()
>>> init_qw = mg.at_node['surface_water__discharge'].copy()
>>> fa = FlowAccumulator(mg)
>>> lmb = LakeMapperBarnes(mg, method='Steepest',
... fill_flat=False, track_lakes=True,
... redirect_flow_steepest_descent=False,
... ignore_overfill=True)
>>> lmb.run_one_step()
>>> np.all(mg.at_node['flow__receiver_node'] == init_flowdirs)
True
>>> lmb = LakeMapperBarnes(mg, method='Steepest',
... fill_flat=False, track_lakes=True,
... redirect_flow_steepest_descent=True,
... ignore_overfill=True)
>>> lmb.run_one_step()
>>> np.all(mg.at_node['flow__receiver_node'] == init_flowdirs)
False
However, note that unless the reaccumulate_flow argument is also set, the ‘drainage_area’ and ‘surface_water__discharge’ fields won’t also get updated:
>>> np.all(mg.at_node['drainage_area'] == init_areas)
True
>>> np.all(mg.at_node['surface_water__discharge'] == init_qw)
True
>>> fa = FlowAccumulator(mg)
>>> lmb = LakeMapperBarnes(mg, method='Steepest',
... fill_flat=False, track_lakes=True,
... redirect_flow_steepest_descent=True,
... reaccumulate_flow=True,
... ignore_overfill=True)
>>> lmb.run_one_step()
>>> np.all(mg.at_node['drainage_area'] == init_areas)
False
>>> np.all(mg.at_node['surface_water__discharge'] == init_qw)
False
Be sure to set both redirect_flow_steepest_descent and reaccumulate_flow to True if you want to reaccumulate flow…
>>> try:
... lmb = LakeMapperBarnes(mg, method='Steepest',
... fill_flat=False, track_lakes=True,
... redirect_flow_steepest_descent=False,
... reaccumulate_flow=True,
... ignore_overfill=True)
... except ValueError:
... print('Oops!')
Oops!
The component is completely happy with irregular grids:
>>> from landlab import HexModelGrid, FieldError
>>> hmg = HexModelGrid((5, 4), spacing=2.)
>>> z_hex = hmg.add_zeros("topographic__elevation", at="node")
>>> z_hex[:] = hmg.node_x
>>> z_hex[11] = -3.
>>> z_hex[12] = -1.
>>> z_hex_init = z_hex.copy()
>>> z_hex
array([ 2., 4., 6., 8.,
1., 3., 5., 7., 9.,
0., 2., -3., -1., 8., 10.,
1., 3., 5., 7., 9.,
2., 4., 6., 8.])
As you can see, nodes 11 and 12 are now a pit. If they were to fill they would fill to the level of 2, the lowest downstream value.
>>> fa = FlowAccumulator(hmg)
>>> lmb = LakeMapperBarnes(
... hmg,
... method='Steepest',
... fill_flat=True,
... track_lakes=False)
>>> lmb.run_one_step()
>>> np.allclose(z_hex[10:13], 2.)
True
>>> hmg = HexModelGrid((5, 4), spacing=2.0)
>>> z_hex = hmg.add_zeros("topographic__elevation", at="node")
>>> z_hex[:] = z_hex_init
>>> try:
... lmb = LakeMapperBarnes(hmg, method='Steepest',
... fill_flat=False,
... surface=z_hex_init,
... redirect_flow_steepest_descent=True,
... track_lakes=True)
... except FieldError:
... print("Oops!") # flowdir field must already exist!
Oops!
>>> fd = FlowDirectorSteepest(hmg)
>>> fa = FlowAccumulator(hmg)
>>> lmb = LakeMapperBarnes(hmg, method='Steepest',
... fill_flat=False, surface=z_hex_init,
... redirect_flow_steepest_descent=True,
... track_lakes=True)
>>> fd.run_one_step()
>>> lmb.run_one_step()
>>> np.allclose(z_hex[10:13], 2.)
True
>>> z_hex[11] > z_hex[10]
True
>>> z_hex[12] > z_hex[11]
True
>>> np.allclose(lmb.lake_depths[10:14], np.array([ 0., 5., 3., 0.]))
True
>>> np.testing.assert_array_almost_equal(
... lmb.lake_volumes,
... 27.712,
... decimal=3)
Together, all this means that we can now run a topographic growth model that permits flooding as it runs:
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import LakeMapperBarnes, FlowAccumulator
>>> from landlab.components import FlowDirectorSteepest
>>> from landlab.components import FastscapeEroder
>>> mg = RasterModelGrid((6, 8))
>>> for edge in ('right', 'top', 'bottom'):
... mg.status_at_node[mg.nodes_at_edge(edge)] = mg.BC_NODE_IS_CLOSED
Because it is what we want the FastscapeEroder to see and work on, it’s actually the water surface that needs to go in as ‘topographic__elevation’. We’ll also need to keep track of the bed elevation though, since the LakeMapper will need it. We start them equal (i.e., topo starts dry).
>>> z_water = mg.add_zeros("topographic__elevation", at="node", dtype=float)
>>> z_water[:] = mg.node_x
>>> z_water[11] = 1.5
>>> z_water[19] = 0.5
>>> z_water[34] = 1.1
>>> z_bed = mg.add_zeros("bedrock__elevation", at="node", dtype=float)
>>> z_bed[:] = z_water # topo starts dry
Let’s just take a look:
>>> np.all(np.equal(
... np.round(z_water, 2),
... np.array([0. , 1. , 2. , 3. , 4. , 5. , 6. , 7. ,
... 0. , 1. , 2. , 1.5, 4. , 5. , 6. , 7. ,
... 0. , 1. , 2. , 0.5, 4. , 5. , 6. , 7. ,
... 0. , 1. , 2. , 3. , 4. , 5. , 6. , 7. ,
... 0. , 1. , 1.1, 3. , 4. , 5. , 6. , 7. ,
... 0. , 1. , 2. , 3. , 4. , 5. , 6. , 7. ])))
True
>>> fd = FlowDirectorSteepest(mg)
>>> fa = FlowAccumulator(mg)
>>> lmb = LakeMapperBarnes(mg, method='D8', fill_flat=True,
... surface='bedrock__elevation',
... fill_surface='topographic__elevation',
... redirect_flow_steepest_descent=True,
... reaccumulate_flow=True,
... track_lakes=True)
>>> sp = FastscapeEroder(mg, K_sp=1., m_sp=0., n_sp=1.)
>>> fd.run_one_step()
>>> fa.run_one_step() # node 18 is draining into the pit...
>>> np.isclose(mg.at_node['topographic__steepest_slope'][18], 1.5)
True
>>> np.allclose(mg.at_node['drainage_area'],
... np.array([ 0., 0., 0., 0., 0., 0., 0., 0.,
... 2., 2., 1., 4., 3., 2., 1., 0.,
... 1., 1., 1., 13., 3., 2., 1., 0.,
... 2., 2., 1., 4., 3., 2., 1., 0.,
... 6., 6., 5., 4., 3., 2., 1., 0.,
... 0., 0., 0., 0., 0., 0., 0., 0.]))
True
>>> lmb.run_one_step() # now node 18 drains correctly, outward ->
>>> np.isclose(mg.at_node['topographic__steepest_slope'][18], 1.)
True
>>> np.allclose(mg.at_node['drainage_area'],
... np.array([ 0., 0., 0., 0., 0., 0., 0., 0.,
... 13., 13., 12., 4., 3., 2., 1., 0.,
... 2., 2., 1., 7., 3., 2., 1., 0.,
... 2., 2., 1., 1., 3., 2., 1., 0.,
... 7., 7., 6., 4., 3., 2., 1., 0.,
... 0., 0., 0., 0., 0., 0., 0., 0.]))
True
>>> np.all(np.equal(
... np.round(z_water, 2),
... np.array([0. , 1. , 2. , 3. , 4. , 5. , 6. , 7. ,
... 0. , 1. , 2. , 2. , 4. , 5. , 6. , 7. ,
... 0. , 1. , 2. , 2. , 4. , 5. , 6. , 7. ,
... 0. , 1. , 2. , 3. , 4. , 5. , 6. , 7. ,
... 0. , 1. , 1.1, 3. , 4. , 5. , 6. , 7. ,
... 0. , 1. , 2. , 3. , 4. , 5. , 6. , 7. ])))
True
>>> sp.run_one_step(0.05) # note m=0 to illustrate effect of slopes
>>> np.all(np.equal(
... np.round(z_water, 2),
... np.array([0. , 1. , 2. , 3. , 4. , 5. , 6. , 7. ,
... 0. , 0.95, 1.95, 2. , 3.9 , 4.95, 5.95, 7. ,
... 0. , 0.95, 1.95, 2. , 3.9 , 4.95, 5.95, 7. ,
... 0. , 0.95, 1.95, 2.93, 3.93, 4.95, 5.95, 7. ,
... 0. , 0.95, 1.09, 2.91, 3.95, 4.95, 5.95, 7. ,
... 0. , 1. , 2. , 3. , 4. , 5. , 6. , 7. ])))
True
If we want to keep this going honouring the depths of the lakes try this next in your loop:
>>> z_bed[:] = np.minimum(z_water, z_bed)
>>> np.all(np.equal(
... np.round(z_bed, 2),
... np.array([0. , 1. , 2. , 3. , 4. , 5. , 6. , 7. ,
... 0. , 0.95, 1.95, 1.5 , 3.9 , 4.95, 5.95, 7. ,
... 0. , 0.95, 1.95, 0.5 , 3.9 , 4.95, 5.95, 7. ,
... 0. , 0.95, 1.95, 2.93, 3.93, 4.95, 5.95, 7. ,
... 0. , 0.95, 1.09, 2.91, 3.95, 4.95, 5.95, 7. ,
... 0. , 1. , 2. , 3. , 4. , 5. , 6. , 7. ])))
True
>>> fd.run_one_step()
>>> fa.run_one_step()
>>> lmb.run_one_step()
Lake node depths are now updated in lmb:
>>> np.round(
... [lmb.lake_depths[lake] for lake in lmb.lake_dict.values()], 2)
array([[ 0.45, 1.45]])
…and the “topography” (i.e., water surface) at the flooded nodes has lowered itself as the lip of the outlet was eroded in the last step:
>>> np.all(np.equal(
... np.round(z_water, 2),
... np.array([0. , 1. , 2. , 3. , 4. , 5. , 6. , 7. ,
... 0. , 0.95, 1.95, 1.95 , 3.9 , 4.95, 5.95, 7. ,
... 0. , 0.95, 1.95, 1.95 , 3.9 , 4.95, 5.95, 7. ,
... 0. , 0.95, 1.95, 2.93, 3.93, 4.95, 5.95, 7. ,
... 0. , 0.95, 1.09, 2.91, 3.95, 4.95, 5.95, 7. ,
... 0. , 1. , 2. , 3. , 4. , 5. , 6. , 7. ])))
True
>>> sp.run_one_step(0.05)
…and so on.
Note that this approach, without passing flooded_nodes to the FastscapeEroder run method, is both more “Landlabbic” and also ensures the information about the lake and the water surface topography are all updated cleanly and correctly.
was_there_overfill¶If the ignore_overfill flag was set to True at instantiation, this property indicates if any depression in the grid has, at any point, been overfilled.
Examples
>>> from landlab import RasterModelGrid
>>> from landlab.components import LakeMapperBarnes, FlowAccumulator
>>> mg = RasterModelGrid((3, 7))
>>> for edge in ('top', 'right', 'bottom'):
... mg.status_at_node[mg.nodes_at_edge(edge)] = mg.BC_NODE_IS_CLOSED
>>> z = mg.add_zeros("topographic__elevation", at="node", dtype=float)
>>> z.reshape(mg.shape)[1, 1:-1] = [1., 0.2, 0.1,
... 1.0000000000000004, 1.5]
>>> z_init = z.copy()
>>> fa = FlowAccumulator(mg)
>>> lmb = LakeMapperBarnes(mg, method='Steepest', fill_flat=True,
... redirect_flow_steepest_descent=False,
... ignore_overfill=False, track_lakes=True)
>>> lmb.run_one_step()
>>> try:
... lmb.was_there_overfill
... except ValueError:
... print('ValueError was raised: ' +
... 'was_there_overfill is only defined if filling to an ' +
... 'inclined surface!')
ValueError was raised: was_there_overfill is only defined if filling to an inclined surface!
>>> z_init = z.copy()
>>> lmb = LakeMapperBarnes(mg, method='Steepest', fill_flat=False,
... redirect_flow_steepest_descent=False,
... ignore_overfill=False, track_lakes=True)
>>> try:
... lmb.run_one_step()
... except ValueError:
... print('ValueError was raised: ' +
... 'Pit is overfilled due to creation of two outlets ' +
... 'as the minimum gradient gets applied. Suppress this ' +
... 'Error with the ignore_overfill flag at component ' +
... 'instantiation.')
ValueError was raised: Pit is overfilled due to creation of two outlets as the minimum gradient gets applied. Suppress this Error with the ignore_overfill flag at component instantiation.
>>> z_init = z.copy()
>>> lmb = LakeMapperBarnes(mg, method='Steepest', fill_flat=False,
... redirect_flow_steepest_descent=False,
... ignore_overfill=True, track_lakes=True)
>>> lmb.run_one_step()
>>> lmb.was_there_overfill
True
>>> z.reshape(mg.shape)[1, 1:-1] = [1., 0.2, 0.1, 1., 1.5]
>>> lmb.run_one_step()
>>> lmb.was_there_overfill # still true as was in the previous example
True
>>> z.reshape(mg.shape)[1, 1:-1] = [1., 0.2, 0.1, 1., 1.5]
>>> lmb = LakeMapperBarnes(mg, method='Steepest', fill_flat=False,
... redirect_flow_steepest_descent=False,
... ignore_overfill=True, track_lakes=True)
>>> lmb.run_one_step()
>>> lmb.was_there_overfill # Now reset
False
>>> lmb.run_one_step() # 2nd run on same fill_surface creates overfill
>>> lmb.was_there_overfill
True
Note however that in this last example, values have NOT changed.
LandslideProbability(grid, number_of_iterations=250, g=9.80665, groundwater__recharge_distribution='uniform', groundwater__recharge_min_value=20.0, groundwater__recharge_max_value=120.0, groundwater__recharge_mean=None, groundwater__recharge_standard_deviation=None, groundwater__recharge_HSD_inputs=[], seed=0)[source]¶Bases: landlab.core.model_component.Component
Landslide probability component using the infinite slope stability model.
Landlab component designed to calculate probability of failure at each grid node based on the infinite slope stability model stability index (Factor of Safety).
The driving force for failure is provided by the user in the form of groundwater recharge; four options for providing recharge are supported. The model uses topographic and soil characteristics provided as input by the user.
The main method of the LandslideProbability class is calculate_landslide_probability()`, which calculates the mean soil relative wetness, probability of soil saturation, and probability of failure at each node based on a Monte Carlo simulation.
Usage:
Option 1 - Uniform recharge
LandslideProbability(grid,
number_of_iterations=250,
groundwater__recharge_distribution='uniform',
groundwater__recharge_min_value=5.,
groundwater__recharge_max_value=121.)
Option 2 - Lognormal recharge
LandslideProbability(grid,
number_of_iterations=250,
groundwater__recharge_distribution='lognormal',
groundwater__recharge_mean=30.,
groundwater__recharge_standard_deviation=0.25)
Option 3 - Lognormal_spatial recharge
LandslideProbability(grid,
number_of_iterations=250,
groundwater__recharge_distribution='lognormal_spatial',
groundwater__recharge_mean=np.random.randint(20, 120, grid_size),
groundwater__recharge_standard_deviation=np.random.rand(grid_size))
Option 4 - Data_driven_spatial recharge
LandslideProbability(grid,
number_of_iterations=250,
groundwater__recharge_distribution='data_driven_spatial',
groundwater__recharge_HSD_inputs=[HSD_dict,
HSD_id_dict,
fract_dict])
Examples
>>> from landlab import RasterModelGrid
>>> from landlab.components.landslides import LandslideProbability
>>> import numpy as np
Create a grid on which to calculate landslide probability.
>>> grid = RasterModelGrid((5, 4), xy_spacing=(0.2, 0.2))
Check the number of core nodes.
>>> grid.number_of_core_nodes
6
The grid will need some input data. To check the names of the fields that provide the input to this component, use the input_var_names class property.
>>> sorted(LandslideProbability.input_var_names) # doctest: +NORMALIZE_WHITESPACE
['soil__density',
'soil__internal_friction_angle',
'soil__maximum_total_cohesion',
'soil__minimum_total_cohesion',
'soil__mode_total_cohesion',
'soil__saturated_hydraulic_conductivity',
'soil__thickness',
'soil__transmissivity',
'topographic__slope',
'topographic__specific_contributing_area']
Check the units for the fields.
>>> LandslideProbability.var_units('topographic__specific_contributing_area')
'm'
Create an input field.
>>> grid.at_node['topographic__slope'] = np.random.rand(grid.number_of_nodes)
If you are not sure about one of the input or output variables, you can get help for specific variables.
>>> LandslideProbability.var_help('soil__transmissivity') # doctest: +NORMALIZE_WHITESPACE
name: soil__transmissivity
description:
mode rate of water transmitted through a unit width of saturated
soil - either provided or calculated with Ksat and soil depth
units: m2/day
unit agnostic: False
at: node
intent: in
Additional required fields for component.
>>> scatter_dat = np.random.randint(1, 10, grid.number_of_nodes)
>>> grid.at_node['topographic__specific_contributing_area'] = np.sort(
... np.random.randint(30, 900, grid.number_of_nodes).astype(float))
>>> grid.at_node['soil__transmissivity'] = np.sort(
... np.random.randint(5, 20, grid.number_of_nodes).astype(float), -1)
>>> grid.at_node['soil__saturated_hydraulic_conductivity'] = np.sort(
... np.random.randint(2, 10, grid.number_of_nodes).astype(float), -1)
>>> grid.at_node['soil__mode_total_cohesion'] = np.sort(
... np.random.randint(30, 900, grid.number_of_nodes).astype(float))
>>> grid.at_node['soil__minimum_total_cohesion'] = (
... grid.at_node['soil__mode_total_cohesion'] - scatter_dat)
>>> grid.at_node['soil__maximum_total_cohesion'] = (
... grid.at_node['soil__mode_total_cohesion'] + scatter_dat)
>>> grid.at_node['soil__internal_friction_angle'] = np.sort(
... np.random.randint(26, 40, grid.number_of_nodes).astype(float))
>>> grid.at_node['soil__thickness'] = np.sort(
... np.random.randint(1, 10, grid.number_of_nodes).astype(float))
>>> grid.at_node['soil__density'] = (2000. * np.ones(grid.number_of_nodes))
Instantiate the ‘LandslideProbability’ component to work on this grid, and run it.
>>> ls_prob = LandslideProbability(grid)
>>> np.allclose(grid.at_node['landslide__probability_of_failure'], 0.)
True
Run the calculate_landslide_probability method to update output variables with grid
>>> ls_prob.calculate_landslide_probability()
Check the output variable names.
>>> sorted(ls_prob.output_var_names) # doctest: +NORMALIZE_WHITESPACE
['landslide__probability_of_failure',
'soil__mean_relative_wetness',
'soil__probability_of_saturation']
Check the output from the component, including array at one node.
>>> np.allclose(grid.at_node['landslide__probability_of_failure'], 0.)
False
>>> core_nodes = ls_prob.grid.core_nodes
References
Required Software Citation(s) Specific to this Component
Strauch, R., Istanbulluoglu, E., Nudurupati, S., Bandaragoda, C., Gasparini, N., Tucker, G. (2018). A hydroclimatological approach to predicting regional landslide probability using Landlab Earth Surface Dynamics 6(1), 49-75. https://dx.doi.org/10.5194/esurf-6-49-2018
Additional References
None Listed
| Parameters: |
|
|---|
__init__(grid, number_of_iterations=250, g=9.80665, groundwater__recharge_distribution='uniform', groundwater__recharge_min_value=20.0, groundwater__recharge_max_value=120.0, groundwater__recharge_mean=None, groundwater__recharge_standard_deviation=None, groundwater__recharge_HSD_inputs=[], seed=0)[source]¶| Parameters: |
|
|---|
calculate_factor_of_safety(i)[source]¶Method to calculate factor of safety.
Method calculates factor-of-safety stability index by using node specific parameters, creating distributions of these parameters, and calculating the index by sampling these distributions ‘n’ times.
The index is calculated from the ‘infinite slope stabilty factor-of-safety equation’ in the format of Pack RT, Tarboton DG, and Goodwin CN (1998),The SINMAP approach to terrain stability mapping.
| Parameters: | i (int) – index of core node ID. |
|---|
calculate_landslide_probability()[source]¶Main method of Landslide Probability class.
Method creates arrays for output variables then loops through all the core nodes to run the method ‘calculate_factor_of_safety.’ Output parameters probability of failure, mean relative wetness, and probability of saturation are assigned as fields to nodes.
LateralEroder(grid, latero_mech='UC', alph=0.8, Kv=0.001, Kl_ratio=1.0, solver='basic', inlet_on=False, inlet_node=None, inlet_area=None, qsinlet=0.0, flow_accumulator=None)[source]¶Bases: landlab.core.model_component.Component
Laterally erode neighbor node through fluvial erosion.
Landlab component that finds a neighbor node to laterally erode and calculates lateral erosion. See the publication:
Langston, A.L., Tucker, G.T.: Developing and exploring a theory for the lateral erosion of bedrock channels for use in landscape evolution models. Earth Surface Dynamics, 6, 1-27, https://doi.org/10.5194/esurf-6-1-2018
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator, LateralEroder
>>> np.random.seed(2010)
Define grid and initial topography
>>> mg = RasterModelGrid((5, 4), xy_spacing=10.0)
>>> mg.set_status_at_node_on_edges(
... right=mg.BC_NODE_IS_CLOSED,
... top=mg.BC_NODE_IS_CLOSED,
... left=mg.BC_NODE_IS_CLOSED,
... bottom=mg.BC_NODE_IS_CLOSED,
... )
>>> mg.status_at_node[1] = mg.BC_NODE_IS_FIXED_VALUE
>>> mg.add_zeros("topographic__elevation", at="node")
array([ 0., 0., 0., 0.,
0., 0., 0., 0.,
0., 0., 0., 0.,
0., 0., 0., 0.,
0., 0., 0., 0.])
>>> rand_noise=np.array(
... [
... 0.00436992, 0.03225985, 0.03107455, 0.00461312,
... 0.03771756, 0.02491226, 0.09613959, 0.07792969,
... 0.08707156, 0.03080568, 0.01242658, 0.08827382,
... 0.04475065, 0.07391732, 0.08221057, 0.02909259,
... 0.03499337, 0.09423741, 0.01883171, 0.09967794,
... ]
... )
>>> mg.at_node["topographic__elevation"] += (
... mg.node_y / 10. + mg.node_x / 10. + rand_noise
... )
>>> U = 0.001
>>> dt = 100
Instantiate flow accumulation and lateral eroder and run each for one step
>>> fa = FlowAccumulator(
... mg,
... surface="topographic__elevation",
... flow_director="FlowDirectorD8",
... runoff_rate=None,
... depression_finder=None,
... )
>>> latero = LateralEroder(mg, latero_mech="UC", Kv=0.001, Kl_ratio=1.5)
Run one step of flow accumulation and lateral erosion to get the dzlat array needed for the next part of the test.
>>> fa.run_one_step()
>>> mg, dzlat = latero.run_one_step(dt)
Evolve the landscape until the first occurence of lateral erosion. Save arrays volume of lateral erosion and topographic elevation before and after the first occurence of lateral erosion
>>> while min(dzlat) == 0.0:
... oldlatvol = mg.at_node["volume__lateral_erosion"].copy()
... oldelev = mg.at_node["topographic__elevation"].copy()
... fa.run_one_step()
... mg, dzlat = latero.run_one_step(dt)
... newlatvol = mg.at_node["volume__lateral_erosion"]
... newelev = mg.at_node["topographic__elevation"]
... mg.at_node["topographic__elevation"][mg.core_nodes] += U * dt
Before lateral erosion occurs, volume__lateral_erosion has values at nodes 6 and 10.
>>> np.around(oldlatvol, decimals=0)
array([ 0., 0., 0., 0.,
0., 0., 79., 0.,
0., 0., 24., 0.,
0., 0., 0., 0.,
0., 0., 0., 0.])
After lateral erosion occurs at node 6, volume__lateral_erosion is reset to 0
>>> np.around(newlatvol, decimals=0)
array([ 0., 0., 0., 0.,
0., 0., 0., 0.,
0., 0., 24., 0.,
0., 0., 0., 0.,
0., 0., 0., 0.])
After lateral erosion at node 6, elevation at node 6 is reduced by -1.41 (the elevation change stored in dzlat[6]). It is also provided as the at-node grid field lateral_erosion__depth_increment.
>>> np.around(oldelev, decimals=2)
array([ 0. , 1.03, 2.03, 3. ,
1.04, 1.77, 2.45, 4.08,
2.09, 2.65, 3.18, 5.09,
3.04, 3.65, 4.07, 6.03,
4.03, 5.09, 6.02, 7.1 ])
>>> np.around(newelev, decimals=2)
array([ 0. , 1.03, 2.03, 3. ,
1.04, 1.77, 1.03, 4.08,
2.09, 2.65, 3.18, 5.09,
3.04, 3.65, 4.07, 6.03,
4.03, 5.09, 6.02, 7.1 ])
>>> np.around(dzlat, decimals=2)
array([ 0. , 0. , 0. , 0. ,
0. , 0. , -1.41, 0. ,
0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. ])
References
Required Software Citation(s) Specific to this Component
Langston, A., Tucker, G. (2018). Developing and exploring a theory for the lateral erosion of bedrock channels for use in landscape evolution models. Earth Surface Dynamics 6(1), 1–27. https://dx.doi.org/10.5194/esurf-6-1-2018
Additional References
None Listed
| Parameters: |
|
|---|
__init__(grid, latero_mech='UC', alph=0.8, Kv=0.001, Kl_ratio=1.0, solver='basic', inlet_on=False, inlet_node=None, inlet_area=None, qsinlet=0.0, flow_accumulator=None)[source]¶| Parameters: |
|
|---|
LinearDiffuser(grid, linear_diffusivity=0.01, method='simple', deposit=True)[source]¶Bases: landlab.core.model_component.Component
This component implements linear diffusion of a Landlab field.
Component assumes grid does not deform. If the boundary conditions on the
grid change after component instantiation, be sure to also call
updated_boundary_conditions to ensure these are reflected in the
component (especially if fixed_links are present).
The method keyword allows control of the way the solver works. Options other than ‘simple’ will make the component run slower, but give second order solutions that incorporate information from more distal nodes (i.e., the diagonals). Note that the option ‘resolve_on_patches’ can result in somewhat counterintuitive behaviour - this option tells the component to treat the diffusivity as a field with directionality to it (i.e., that the diffusivites are defined on links). Thus if all links have the same diffusivity value, with this flag active “effective” diffusivities at the nodes will be higher than this value (by a factor of root 2) as the diffusivity at each patch will be the mean vector sum of that at the bounding links.
The primary method of this class is run_one_step.
Examples
>>> from landlab import RasterModelGrid
>>> import numpy as np
>>> mg = RasterModelGrid((9, 9))
>>> z = mg.add_zeros("topographic__elevation", at="node")
>>> z.reshape((9, 9))[4, 4] = 1.
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, True)
>>> ld = LinearDiffuser(mg, linear_diffusivity=1.)
>>> for i in range(1):
... ld.run_one_step(1.)
>>> np.isclose(z[mg.core_nodes].sum(), 1.)
True
>>> mg2 = RasterModelGrid((5, 30))
>>> z2 = mg2.add_zeros("topographic__elevation", at="node")
>>> z2.reshape((5, 30))[2, 8] = 1.
>>> z2.reshape((5, 30))[2, 22] = 1.
>>> mg2.set_closed_boundaries_at_grid_edges(True, True, True, True)
>>> kd = mg2.node_x/mg2.node_x.mean()
>>> ld2 = LinearDiffuser(mg2, linear_diffusivity=kd)
>>> for i in range(10):
... ld2.run_one_step(0.1)
>>> z2[mg2.core_nodes].sum() == 2.
True
>>> z2.reshape((5, 30))[2, 8] > z2.reshape((5, 30))[2, 22]
True
An example using links:
>>> mg1 = RasterModelGrid((10, 10), xy_spacing=100.)
>>> mg2 = RasterModelGrid((10, 10), xy_spacing=100.)
>>> z1 = mg1.add_zeros("topographic__elevation", at="node")
>>> z2 = mg2.add_zeros("topographic__elevation", at="node")
>>> dt = 1.
>>> nt = 10
>>> kappa_links = mg2.add_ones("surface_water__discharge", at="link")
>>> kappa_links *= 10000.
>>> dfn1 = LinearDiffuser(mg1, linear_diffusivity=10000.)
>>> dfn2 = LinearDiffuser(mg2, linear_diffusivity='surface_water__discharge')
>>> for i in range(nt):
... z1[mg1.core_nodes] += 1.
... z2[mg2.core_nodes] += 1.
... dfn1.run_one_step(dt)
... dfn2.run_one_step(dt)
>>> np.allclose(z1, z2)
True
>>> z2.fill(0.)
>>> dfn2 = LinearDiffuser(mg2, linear_diffusivity='surface_water__discharge',
... method='resolve_on_patches')
>>> for i in range(nt):
... z2[mg2.core_nodes] += 1.
... dfn2.run_one_step(dt)
>>> np.all(z2[mg2.core_nodes] < z1[mg2.core_nodes])
True
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
Culling, W. (1963). Soil Creep and the Development of Hillside Slopes. The Journal of Geology 71(2), 127-161. https://dx.doi.org/10.1086/626891
| Parameters: |
|
|---|
__init__(grid, linear_diffusivity=0.01, method='simple', deposit=True)[source]¶| Parameters: |
|
|---|
fixed_grad_anchors¶Fixed gradient anchors.
fixed_grad_nodes¶Fixed gradient nodes.
fixed_grad_offsets¶Fixed gradient offsets.
run_one_step(dt)[source]¶Run the diffuser for one timestep, dt.
If the imposed timestep dt is longer than the Courant-Friedrichs-Lewy condition for the diffusion, this timestep will be internally divided as the component runs, as needed.
| Parameters: | dt (float (time)) – The imposed timestep. |
|---|
time_step¶Returns internal time-step size (as a property).
updated_boundary_conditions()[source]¶Call if grid BCs are updated after component instantiation.
Sets fixed_grad_nodes, fixed_grad_anchors, & fixed_grad_offsets, such that:
value[fixed_grad_nodes] = value[fixed_grad_anchors] + offset
Examples
>>> from landlab import RasterModelGrid
>>> import numpy as np
>>> mg = RasterModelGrid((4, 5))
>>> z = mg.add_zeros("topographic__elevation", at="node")
>>> z[mg.core_nodes] = 1.
>>> ld = LinearDiffuser(mg, linear_diffusivity=1.)
>>> ld.fixed_grad_nodes.size == 0
True
>>> ld.fixed_grad_anchors.size == 0
True
>>> ld.fixed_grad_offsets.size == 0
True
>>> mg.at_link['topographic__slope'] = mg.calc_grad_at_link(
... 'topographic__elevation')
>>> mg.status_at_node[mg.perimeter_nodes] = mg.BC_NODE_IS_FIXED_GRADIENT
>>> ld.updated_boundary_conditions()
>>> ld.fixed_grad_nodes
array([ 1, 2, 3, 5, 9, 10, 14, 16, 17, 18])
>>> ld.fixed_grad_anchors
array([ 6, 7, 8, 6, 8, 11, 13, 11, 12, 13])
>>> ld.fixed_grad_offsets
array([-1., -1., -1., -1., -1., -1., -1., -1., -1., -1.])
>>> np.allclose(z[ld.fixed_grad_nodes],
... z[ld.fixed_grad_anchors] + ld.fixed_grad_offsets)
True
ListricKinematicExtender(grid, extension_rate=0.001, fault_dip=60.0, fault_location=None, detachment_depth=10000.0, track_crustal_thickness=False, fields_to_shift=[])[source]¶Bases: landlab.core.model_component.Component
Apply tectonic extension and subsidence kinematically to a raster or hex grid.
Extension is east-west, with a north-south fault in the case of a raster grid, and a N30E fault in the case of a hex grid (obique extension).
Examples
>>> from landlab import RasterModelGrid
>>> from landlab.components import ListricKinematicExtender
>>> grid = RasterModelGrid((3, 7), xy_spacing=2500.0)
>>> topo = grid.add_zeros('topographic__elevation', at='node')
>>> lke = ListricKinematicExtender(grid, fault_location=2500.0)
>>> lke.update_subsidence_rate()
>>> np.round(grid.at_node["subsidence_rate"][8:13], 6)
array([ 0. , 0.001123, 0.000729, 0.000472, 0.000306])
Deform vertically and horizontally to represent tectonic extension.
| Parameters: |
|
|---|
__init__(grid, extension_rate=0.001, fault_dip=60.0, fault_location=None, detachment_depth=10000.0, track_crustal_thickness=False, fields_to_shift=[])[source]¶Deform vertically and horizontally to represent tectonic extension.
| Parameters: |
|
|---|
fault_plane_elev(dist_from_fault_trace)[source]¶Return elevation of fault plane at given distance from fault trace.
Examples
>>> from landlab import RasterModelGrid
>>> from landlab.components import ListricKinematicExtender
>>> grid = RasterModelGrid((3, 7), xy_spacing=5000.0)
>>> topo = grid.add_zeros('topographic__elevation', at='node')
>>> lke = ListricKinematicExtender(grid, fault_location=10000.0)
>>> np.round(lke.fault_plane_elev(np.array([0.0, 5000.0, 10000.0])))
array([ -0., -5794., -8231.])
LithoLayers(grid, z0s, ids, attrs, x0=0, y0=0, function=<function LithoLayers.<lambda>>, layer_type='EventLayers', dz_advection=0, rock_id=None)[source]¶Bases: landlab.components.lithology.lithology.Lithology
Create LithoLayers component.
A LithoLayers is a three dimentional representation of material operated on by landlab components. Material can be removed through erosion or added to through deposition. Rock types can have multiple attributes (e.g. age, erodability or other parameter values, etc).
If the tracked properties are model grid fields, they will be updated to the surface values of the Lithology. If the properties are not grid fields then at-node grid fields will be created with their names.
It is constructed by specifying a series of depths below the surface, an anchor point, a series of rock type ids, and the functional form of a surface. Depths and IDs are both specified in order of closest to the surface to furthest from the surface.
Additionally, an attribute dictionary specifies the properties of each rock type. This dictionary is expected to have the form of:
attrs = {'K_sp': {1: 0.001,
2: 0.0001},
'D': {1: 0.01,
2: 0.001}}
Where 'K_sp' and 'D' are properties to track, and 1 and 2
are rock type IDs. The rock type IDs can be any type that is valid as a
python dictionary key.
References
Required Software Citation(s) Specific to this Component
Barnhart, K., Hutton, E., Gasparini, N., Tucker, G. (2018). Lithology: A Landlab submodule for spatially variable rock properties. Journal of Open Source Software 3(30), 979 - 2. https://dx.doi.org/10.21105/joss.00979
Additional References
None Listed
Create a new instance of a LithoLayers.
| Parameters: |
|
|---|