landlab.components.flow_accum.lossy_flow_accumulator

Accumulate flow and calc drainage area, while permitting gain or loss of discharge during flow.

DEJH, late 2018

class LossyFlowAccumulator[source]

Bases: FlowAccumulator

Component to calculate drainage area and accumulate flow, while permitting dynamic loss or gain of flow downstream.

This component is closely related to the FlowAccumulator, in that this is accomplished by first finding flow directions by a user-specified method and then calculating the drainage area and discharge. However, this component additionally requires the passing of a function that describes how discharge is lost or gained downstream:

f(Qw, nodeID, linkID, grid)

See the Examples below to see how this works in practice.

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 that the DepressionFinderAndRouter is not particularly intelligent when running on lossy streams, and in particular, it will reroute flow around pits even when they are in fact not filled due to loss.

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.

LossyFlowAccumulator stores as ModelGrid fields:

  • Node array of drainage areas: "drainage_area"

  • Node array of discharges: "surface_water__discharge"

  • Node array of discharge loss in transit (vol / sec). This is the total loss across all of the downstream links: "surface_water__discharge_loss"

  • 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; see the FlowAccumulator for full details. These are:

  • Node array of receivers (nodes that receive flow), or ITS OWN ID if there is no receiver: "flow__receiver_node".

  • Node array of flow proportions: "flow__receiver_proportions".

  • Node array of links carrying flow: "flow__link_to_receiver_node".

  • Node array of downhill slopes from each receiver: "topographic__steepest_slope".

  • Boolean node array of all local lows: "flow__sink_flag".

The primary method of this class is run_one_step.

Examples

These examples pertain only to the LossyFlowAccumulator. See the main FlowAccumulator documentation for more generic and comprehensive examples.

First, a very simple example. Here’s a 50% loss of discharge every time flow moves along a node:

>>> import numpy as np
>>> from landlab import RasterModelGrid, HexModelGrid
>>> from landlab.components import FlowDirectorSteepest
>>> from landlab.components import DepressionFinderAndRouter

>>> mg = RasterModelGrid((3, 5), xy_spacing=(2, 1))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, False, True)
>>> z = mg.add_field("topographic__elevation", mg.node_x + mg.node_y, at="node")

>>> def mylossfunction(qw):
...     return 0.5 * qw
...

>>> fa = LossyFlowAccumulator(
...     mg,
...     "topographic__elevation",
...     flow_director=FlowDirectorSteepest,
...     loss_function=mylossfunction,
... )
>>> fa.run_one_step()

>>> mg.at_node["drainage_area"].reshape(mg.shape)
array([[0., 0., 0., 0., 0.],
       [6., 6., 4., 2., 0.],
       [0., 0., 0., 0., 0.]])
>>> mg.at_node["surface_water__discharge"].reshape(mg.shape)
array([[0.  , 0.  , 0.  , 0.  , 0.  ],
       [1.75, 3.5 , 3.  , 2.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  ]])
>>> mg.at_node["surface_water__discharge_loss"].reshape(mg.shape)
array([[0.  , 0.  , 0.  , 0.  , 0.  ],
       [0.  , 1.75, 1.5 , 1.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  ]])

Here we use a spatially distributed field to derive loss terms, and also use a filled, non-raster grid.

>>> dx = (2.0 / (3.0**0.5)) ** 0.5  # area to be 100.0
>>> hmg = HexModelGrid((5, 3), spacing=dx, xy_of_lower_left=(-1.0745, 0.0))
>>> z = hmg.add_field(
...     "topographic__elevation",
...     hmg.node_x**2 + np.round(hmg.node_y) ** 2,
...     at="node",
... )
>>> z[9] = -10.0  # poke a hole
>>> lossy = hmg.add_zeros("mylossterm", dtype=float, at="node")
>>> lossy[14] = 1.0  # suppress all flow from node 14

Without loss looks like this:

>>> fa = LossyFlowAccumulator(
...     hmg,
...     "topographic__elevation",
...     flow_director=FlowDirectorSteepest,
...     depression_finder=DepressionFinderAndRouter,
... )
>>> fa.run_one_step()
>>> hmg.at_node["flow__receiver_node"]
array([ 0,  1,  2,
        3,  0,  9,  6,
        7,  9,  4,  9, 11,
       12,  9,  9, 15,
       16, 17, 18])
>>> np.round(hmg.at_node["drainage_area"])
array([7., 0., 0.,
       0., 7., 1., 0.,
       0., 1., 6., 1., 0.,
       0., 1., 1., 0.,
       0., 0., 0.])
>>> np.round(hmg.at_node["surface_water__discharge"])
array([7., 0., 0.,
       0., 7., 1., 0.,
       0., 1., 6., 1., 0.,
       0., 1., 1., 0.,
       0., 0., 0.])

With loss looks like this:

>>> def mylossfunction2(Qw, nodeID, linkID, grid):
...     return (1.0 - grid.at_node["mylossterm"][nodeID]) * Qw
...
>>> fa = LossyFlowAccumulator(
...     hmg,
...     "topographic__elevation",
...     flow_director=FlowDirectorSteepest,
...     depression_finder=DepressionFinderAndRouter,
...     loss_function=mylossfunction2,
... )
>>> fa.run_one_step()
>>> np.round(hmg.at_node["drainage_area"])
array([7., 0., 0.,
       0., 7., 1., 0.,
       0., 1., 6., 1., 0.,
       0., 1., 1., 0.,
       0., 0., 0.])
>>> np.round(hmg.at_node["surface_water__discharge"])
array([6., 0., 0.,
       0., 6., 1., 0.,
       0., 1., 5., 1., 0.,
       0., 1., 1., 0.,
       0., 0., 0.])
>>> np.allclose(
...     hmg.at_node["surface_water__discharge_loss"],
...     lossy * hmg.at_node["surface_water__discharge"],
... )
True

(Loss is only happening from the node, 14, that we set it to happen at.)

Finally, note we can use the linkIDs to create flow-length-dependent effects:

>>> from landlab.components import FlowDirectorMFD
>>> mg = RasterModelGrid((4, 6), xy_spacing=(1, 2))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, False, True)
>>> z = mg.add_field("topographic__elevation", 2.0 * mg.node_x, at="node")
>>> z[9] = 8.0
>>> z[16] = 6.5  # force the first node sideways

>>> L = mg.add_zeros("spatialloss", at="node")
>>> mg.at_node["spatialloss"][9] = 1.0
>>> mg.at_node["spatialloss"][13] = 1.0
>>> def fancyloss(Qw, nodeID, linkID, grid):
...     # now a true transmission loss:
...     Lt = 1.0 - 1.0 / grid.length_of_link[linkID] ** 2
...     Lsp = grid.at_node["spatialloss"][nodeID]
...     return Qw * (1.0 - Lt) * (1.0 - Lsp)
...

>>> fa = LossyFlowAccumulator(
...     mg,
...     "topographic__elevation",
...     flow_director=FlowDirectorMFD,
...     loss_function=fancyloss,
... )
>>> fa.run_one_step()

>>> mg.at_node["drainage_area"].reshape(mg.shape)
array([[ 0. ,  0. , 0. ,  0. ,  0. ,  0. ],
       [ 5.6,  5.6, 3.6,  2. ,  2. ,  0. ],
       [10.4, 10.4, 8.4,  6.4,  4. ,  0. ],
       [ 0. ,  0. , 0. ,  0. ,  0. ,  0. ]])
>>> mg.at_node["surface_water__discharge"].reshape(mg.shape)
array([[0. , 0. , 0. , 0. , 0. , 0. ],
       [4. , 4. , 2. , 2. , 2. , 0. ],
       [0. , 8.5, 6.5, 4.5, 2.5, 0. ],
       [0. , 0. , 0. , 0. , 0. , 0. ]])

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.

Parameters:

  • grid (ModelGrid) – A Landlab grid.

  • surface (str, int or array of int) – The surface to direct flow across. Can be a field name at node or an array of length node.

  • flow_director (str, class or instance of class.) – A string of method or class name (e.g. 'D8' or 'FlowDirectorD8'), an uninstantiated FlowDirector class, or an instance of a FlowDirector class. This sets the method used to calculate flow directions.

  • runoff_rate (field name, array, or float, optional (m/time)) – If provided, sets the runoff rate and will be assigned to the grid field 'water__unit_flux_in'. If a spatially and and temporally variable runoff rate is desired, pass this field name and update the field through model run time. If both the field and argument are present at the time of initialization, runoff_rate will overwrite the field. If neither are set, defaults to spatially constant unit input.

  • depression_finder (str, class, instance of class, optional) – A string of class name (e.g., 'DepressionFinderAndRouter'), an uninstantiated DepressionFinder class, or an instance of a DepressionFinder class. This sets the method for depression finding.

  • loss_function (func, optional) –

    A function of the form f(Qw, [node_ID, [linkID, [grid]]]), where Qw is the discharge at a node, node_ID the ID of the node at which the loss is to be calculated, linkID is the ID of the link down which the outflow drains (or a d8 ID if the routing is d8), and grid is a Landlab ModelGrid. The function then returns the new discharge at the node after the function is applied.

    Note that if a linkID is needed, a nodeID must also be specified, even if only as a dummy parameter; similarly, if a grid is to be passed, all of the preceding parameters must be specified. Both nodeID and linkID are required to permit spatially variable losses, and also losses dependent on flow path geometry (e.g., flow length). The grid is passed to allow fields or grid properties describing values across the grid to be accessed for the loss calculation (see examples). This function expects (float, [int, [int, [ModelGrid]]]), and return a single float, the new discharge value. This behavior is verified during component instantiation.

  • **kwargs (optional) – Any additional parameters to pass to a FlowDirector or DepressionFinderAndRouter instance (e.g., partion_method for FlowDirectorMFD). This will have no effect if an instantiated component is passed using the flow_director or depression_finder keywords.

__init__(grid, surface='topographic__elevation', flow_director='FlowDirectorSteepest', runoff_rate=None, depression_finder=None, loss_function=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.

Parameters:

  • grid (ModelGrid) – A Landlab grid.

  • surface (str, int or array of int) – The surface to direct flow across. Can be a field name at node or an array of length node.

  • flow_director (str, class or instance of class.) – A string of method or class name (e.g. 'D8' or 'FlowDirectorD8'), an uninstantiated FlowDirector class, or an instance of a FlowDirector class. This sets the method used to calculate flow directions.

  • runoff_rate (field name, array, or float, optional (m/time)) – If provided, sets the runoff rate and will be assigned to the grid field 'water__unit_flux_in'. If a spatially and and temporally variable runoff rate is desired, pass this field name and update the field through model run time. If both the field and argument are present at the time of initialization, runoff_rate will overwrite the field. If neither are set, defaults to spatially constant unit input.

  • depression_finder (str, class, instance of class, optional) – A string of class name (e.g., 'DepressionFinderAndRouter'), an uninstantiated DepressionFinder class, or an instance of a DepressionFinder class. This sets the method for depression finding.

  • loss_function (func, optional) –

    A function of the form f(Qw, [node_ID, [linkID, [grid]]]), where Qw is the discharge at a node, node_ID the ID of the node at which the loss is to be calculated, linkID is the ID of the link down which the outflow drains (or a d8 ID if the routing is d8), and grid is a Landlab ModelGrid. The function then returns the new discharge at the node after the function is applied.

    Note that if a linkID is needed, a nodeID must also be specified, even if only as a dummy parameter; similarly, if a grid is to be passed, all of the preceding parameters must be specified. Both nodeID and linkID are required to permit spatially variable losses, and also losses dependent on flow path geometry (e.g., flow length). The grid is passed to allow fields or grid properties describing values across the grid to be accessed for the loss calculation (see examples). This function expects (float, [int, [int, [ModelGrid]]]), and return a single float, the new discharge value. This behavior is verified during component instantiation.

  • **kwargs (optional) – Any additional parameters to pass to a FlowDirector or DepressionFinderAndRouter instance (e.g., partion_method for FlowDirectorMFD). This will have no effect if an instantiated component is passed using the flow_director or depression_finder keywords.

static __new__(cls, *args, **kwds)
accumulate_flow(update_flow_director=True, update_depression_finder=True)

Function to make FlowAccumulator calculate drainage area and discharge.

Running accumulate_flow() results in the following to occur:

  1. Flow directions are updated (unless update_flow_director is set as False). This incluldes checking for updated boundary conditions.

  2. The depression finder, if present is updated (unless update_depression_finder is set as False).

  3. Intermediate steps that analyse the drainage network topology and create datastructures for efficient drainage area and discharge calculations.

  4. Calculation of drainage area and discharge.

  5. Return of drainage area and discharge.

Parameters:

  • update_flow_director (optional, bool) – Whether to update the flow director. Default is True.

  • update_depression_finder (optional, bool) – Whether to update the depression finder, if present. Default is True.

Returns:

  • drainage_area (array) – At node array which points to the field grid.at_node[“drainage_area”].

  • surface_water__discharge – At node array which points to the field grid.at_node[“surface_water__discharge”].

cite_as = ''
property coords

Return the coordinates of nodes on grid attached to the component.

property current_time

Current time.

Some components may keep track of the current time. In this case, the current_time attribute is incremented. Otherwise it is set to None.

Return type:

current_time

definitions = (('drainage_area', "Upstream accumulated surface area contributing to the node's discharge"), ('flow__data_structure_delta', "Node array containing the elements delta[1:] of the data structure 'delta' used for construction of the downstream-to-upstream node array"), ('flow__upstream_node_order', 'Node array containing downstream-to-upstream ordered list of node IDs'), ('surface_water__discharge', 'Volumetric discharge of surface water'), ('surface_water__discharge_loss', 'Total volume of water per second lost during all flow out of the node'), ('topographic__elevation', 'Land surface topographic elevation'), ('water__unit_flux_in', 'External volume water per area per time input to each node (e.g., rainfall rate)'))
property depression_finder

The DepressionFinder used internally.

depression_handler_raster_direction_method()

Return ‘D8’ or ‘D4’ depending on the direction method used.

(Note: only call this function for a raster gird; does not handle multiple-flow directors)

flooded_nodes_present()
property flow_director

The FlowDirector used internally.

flow_director_raster_method()

Return ‘D8’ or ‘D4’ depending on the direction method used.

(Note: only call this function for a raster gird; does not handle multiple-flow directors)

classmethod from_path(grid, path)

Create a component from an input file.

Parameters:

  • grid (ModelGrid) – A landlab grid.

  • path (str or file_like) – Path to a parameter file, contents of a parameter file, or a file-like object.

Returns:

A newly-created component.

Return type:

Component

property grid

Return the grid attached to the component.

headwater_nodes()

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]))
initialize_optional_output_fields()

Create fields for a component based on its optional field outputs, if declared in _optional_var_names.

This method will create new fields (without overwrite) for any fields output by the component as optional. New fields are initialized to zero. New fields are created as arrays of floats, unless the component also contains the specifying property _var_type.

initialize_output_fields(values_per_element=None)

Create fields for a component based on its input and output var names.

This method will create new fields (without overwrite) for any fields output by, but not supplied to, the component. New fields are initialized to zero. Ignores optional fields. New fields are created as arrays of floats, unless the component specifies the variable type.

Parameters:

values_per_element (int (optional)) – On occasion, it is necessary to create a field that is of size (n_grid_elements, values_per_element) instead of the default size (n_grid_elements,). Use this keyword argument to acomplish this task.

input_var_names = ()

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)
>>> np.flipud(
...     mg.add_field(
...         "topographic__elevation",
...         mg.node_x + mg.node_y,
...         at="node",
...     ).reshape(mg.shape)
... )
array([[4., 5., 6., 7., 8.],
       [3., 4., 5., 6., 7.],
       [2., 3., 4., 5., 6.],
       [1., 2., 3., 4., 5.],
       [0., 1., 2., 3., 4.]])
>>> fa = FlowAccumulator(mg, "topographic__elevation", flow_director="MFD")
>>> fa.run_one_step()
>>> link_order = fa.link_order_upstream()
>>> link_order  
array([ 5, 14, 10,  6, 11,  7, 23, 19, 15, 20, 16, 28, 24, 29, 25])
>>> link_order[0]
5
>>> sorted(link_order[1:4])
[6, 10, 14]
>>> sorted(link_order[4:9])
[7, 11, 15, 19, 23]
>>> sorted(link_order[9:13])
[16, 20, 24, 28]
>>> sorted(link_order[13:])
[25, 29]
>>> np.all(sorted(link_order) == mg.active_links)
True
name = 'LossyFlowAccumulator'
property node_drainage_area

Return the drainage area.

property node_order_upstream

Return the upstream node order (drainage stack).

property node_water_discharge

Return the surface water discharge.

optional_var_names = ('topographic__elevation', 'water__unit_flux_in')
output_var_names = ('drainage_area', 'flow__data_structure_delta', 'flow__upstream_node_order', 'surface_water__discharge', 'surface_water__discharge_loss')
pits_present()
run_one_step()

Accumulate flow and save to the model grid.

  1. Flow directions are updated. This incluldes checking for updated boundary conditions.

  2. The depression finder, if present is updated.

  3. Intermediate steps that analyse the drainage network topology and create datastructures for efficient drainage area and discharge calculations.

  4. Calculation of drainage area and discharge.

  5. Return of drainage area and discharge.

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.

property shape

Return the grid shape attached to the component, if defined.

property surface_values

Values of the surface over which flow is accumulated.

unit_agnostic = True
units = (('drainage_area', 'm**2'), ('flow__data_structure_delta', '-'), ('flow__upstream_node_order', '-'), ('surface_water__discharge', 'm**3/s'), ('surface_water__discharge_loss', 'm**3/s'), ('topographic__elevation', 'm'), ('water__unit_flux_in', 'm/s'))
classmethod var_definition(name)

Get a description of a particular field.

Parameters:

name (str) – A field name.

Returns:

A description of each field.

Return type:

tuple of (name, *description*)

classmethod var_help(name)

Print a help message for a particular field.

Parameters:

name (str) – A field name.

classmethod var_loc(name)

Location where a particular variable is defined.

Parameters:

name (str) – A field name.

Returns:

The location (‘node’, ‘link’, etc.) where a variable is defined.

Return type:

str

var_mapping = (('drainage_area', 'node'), ('flow__data_structure_delta', 'node'), ('flow__upstream_node_order', 'node'), ('surface_water__discharge', 'node'), ('surface_water__discharge_loss', 'node'), ('topographic__elevation', 'node'), ('water__unit_flux_in', 'node'))
classmethod var_type(name)

Returns the dtype of a field (float, int, bool, str…).

Parameters:

name (str) – A field name.

Returns:

The dtype of the field.

Return type:

dtype

classmethod var_units(name)

Get the units of a particular field.

Parameters:

name (str) – A field name.

Returns:

Units for the given field.

Return type:

str