LakeMapperBarnes: Component to temporarily fill depressions and reroute flow across them#

class LakeMapperBarnes(*args, **kwds)[source]#

Bases: 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:

  • grid (ModelGrid) – A grid.

  • surface (field name at node or array of length node) – The surface to direct flow across.

  • method ({'Steepest', 'D8'}) – Whether or not to recognise diagonals as valid flow paths, if a raster. Otherwise, no effect.

  • fill_flat (bool) – If True, pits will be filled to perfectly horizontal. If False, the new surface will be slightly inclined to give steepest descent flow paths to the outlet.

  • fill_surface (bool) – Sets the field or array to fill. If fill_surface is surface, this operation occurs in place, and is faster. Note that the component will overwrite fill_surface if it exists; to supply an existing water level to it, supply that water level field as surface, not fill_surface.

  • redirect_flow_steepest_descent (bool) – If True, the component outputs modified versions of the ‘flow__receiver_node’, ‘flow__link_to_receiver_node’, ‘flow__sink_flag’, and ‘topographic__steepest_slope’ fields. These are the fields output by the FlowDirector components, so set to True if you wish to pass this LakeFiller to the FlowAccumulator, or if you wish to work directly with the new, correct flow directions and slopes without rerunning these components on your new surface. Ensure the necessary fields already exist, and have already been calculated by a FlowDirector! This also means you need to instantiate your FlowDirector before you instantiate the LakeMapperBarnes. Note that the new topographic__steepest_slope will always be set to zero, even if fill_flat=False (i.e., there is actually a miniscule gradient on the new topography, which gets ignored).

  • reaccumulate_flow (bool) – If True, and redirect_flow_steepest_descent is True, the run method will (re-)accumulate the flow after redirecting the flow. This means the ‘drainage_area’, ‘surface_water__discharge’, ‘flow__upstream_node_order’, and the other various flow accumulation fields (see output field names) will now reflect the new drainage patterns without having to manually reaccumulate the discharge. If True but redirect_flow_steepest_descent is False, raises an ValueError.

  • ignore_overfill (bool) – If True, suppresses the Error that would normally be raised during creation of a gentle incline on a fill surface (i.e., if not fill_flat). Typically this would happen on a synthetic DEM where more than one outlet is possible at the same elevation. If True, the was_there_overfill property can still be used to see if this has occurred.

  • track_lakes (bool) – If True, the component permits a slight hit to performance in order to explicitly track which nodes have been filled, and to enable queries on that data in retrospect. Set to False to simply fill the surface and be done with it.

__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:

  • grid (ModelGrid) – A grid.

  • surface (field name at node or array of length node) – The surface to direct flow across.

  • method ({'Steepest', 'D8'}) – Whether or not to recognise diagonals as valid flow paths, if a raster. Otherwise, no effect.

  • fill_flat (bool) – If True, pits will be filled to perfectly horizontal. If False, the new surface will be slightly inclined to give steepest descent flow paths to the outlet.

  • fill_surface (bool) – Sets the field or array to fill. If fill_surface is surface, this operation occurs in place, and is faster. Note that the component will overwrite fill_surface if it exists; to supply an existing water level to it, supply that water level field as surface, not fill_surface.

  • redirect_flow_steepest_descent (bool) – If True, the component outputs modified versions of the ‘flow__receiver_node’, ‘flow__link_to_receiver_node’, ‘flow__sink_flag’, and ‘topographic__steepest_slope’ fields. These are the fields output by the FlowDirector components, so set to True if you wish to pass this LakeFiller to the FlowAccumulator, or if you wish to work directly with the new, correct flow directions and slopes without rerunning these components on your new surface. Ensure the necessary fields already exist, and have already been calculated by a FlowDirector! This also means you need to instantiate your FlowDirector before you instantiate the LakeMapperBarnes. Note that the new topographic__steepest_slope will always be set to zero, even if fill_flat=False (i.e., there is actually a miniscule gradient on the new topography, which gets ignored).

  • reaccumulate_flow (bool) – If True, and redirect_flow_steepest_descent is True, the run method will (re-)accumulate the flow after redirecting the flow. This means the ‘drainage_area’, ‘surface_water__discharge’, ‘flow__upstream_node_order’, and the other various flow accumulation fields (see output field names) will now reflect the new drainage patterns without having to manually reaccumulate the discharge. If True but redirect_flow_steepest_descent is False, raises an ValueError.

  • ignore_overfill (bool) – If True, suppresses the Error that would normally be raised during creation of a gentle incline on a fill surface (i.e., if not fill_flat). Typically this would happen on a synthetic DEM where more than one outlet is possible at the same elevation. If True, the was_there_overfill property can still be used to see if this has occurred.

  • track_lakes (bool) – If True, the component permits a slight hit to performance in order to explicitly track which nodes have been filled, and to enable queries on that data in retrospect. Set to False to simply fill the surface and be done with it.

property 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.0  # 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.0, 1.0]))
True
property 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.0  # 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(
...     [
...         [0, 0, 0, 0, 0, 0],
...         [0, 1, 0, 1, 0, 0],
...         [0, 0, 1, 1, 0, 0],
...         [0, 0, 0, 0, 1, 0],
...         [0, 0, 0, 0, 0, 0],
...     ],
...     dtype=bool,
... )
>>> np.all(lmb.lake_at_node.reshape(mg.shape) == lake_at_node)
True
property 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.0  # 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()

This won’t work as surface & fill_surface are both z

>>> lmb.lake_depths
Traceback (most recent call last):
...
ValueError: surface and fill_surface must be different fields or arrays
to enable the property lake_depths!

ValueError was raised because 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 = [
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
...     [0.0, 1.0, 0.0, 0.5, 0.0, 0.0],
...     [0.0, 0.0, 0.4, 0.7, 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],
... ]
>>> np.all(
...     np.equal(lmb.lake_depths.reshape(mg.shape), lake_depths)
... )  # slope applied, so...
False
>>> np.allclose(lmb.lake_depths.reshape(mg.shape), lake_depths)
True
property 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.0, 1.0, 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]
>>> mg.at_node["topographic__elevation"] = [
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
...     [0.0, 2.1, 1.1, 0.6, 1.6, 0.0],
...     [0.0, 2.0, 1.0, 0.5, 1.5, 0.0],
...     [0.0, 2.2, 1.2, 0.7, 1.7, 0.0],
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
... ]
>>> z = mg.at_node["topographic__elevation"]
>>> z_init = mg.at_node["topographic__elevation"].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
property 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.0  # 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 = [
...     [-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(lmb.lake_map.reshape(mg.shape) == lake_map)
True

Note that the outlet node is NOT part of the lake.

Updating the elevations works fine:

>>> z.reshape(mg.shape)[1:4, 1:-1] = [
...     [2.1, 1.1, 0.6, 1.6],
...     [2.0, 1.0, 0.5, 1.5],
...     [2.2, 1.2, 0.7, 1.7],
... ]

>>> lmb.run_one_step()
>>> new_lake_map = [
...     [-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(lmb.lake_map.reshape(mg.shape) == new_lake_map)
True
property 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
...
>>> mg.at_node["topographic__elevation"] = [
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
...     [0.0, 2.1, 1.1, 0.6, 1.6, 0.0],
...     [0.0, 2.0, 1.0, 0.5, 1.5, 0.0],
...     [0.0, 2.2, 1.2, 0.7, 1.7, 0.0],
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
... ]
>>> z = mg.at_node["topographic__elevation"]
>>> z_init = mg.at_node["topographic__elevation"].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
property 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.0  # 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()
>>> lmb.lake_volumes  # won't work as surface & fill_surface are both z
Traceback (most recent call last):
...
ValueError: 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()
>>> lmb.lake_outlets
[16, 8]
>>> np.allclose(lmb.lake_volumes, np.array([1.7, 1.0]))
True
property 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.0  # 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.0)
>>> for edge in ("left", "top", "bottom"):
...     mg.status_at_node[mg.nodes_at_edge(edge)] = mg.BC_NODE_IS_CLOSED
...
>>> mg.at_node["topographic__elevation"] = [
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
...     [0.0, 2.1, 1.1, 0.6, 1.6, 0.0],
...     [0.0, 2.0, 1.0, 0.5, 1.5, 0.0],
...     [0.0, 2.2, 1.2, 0.7, 1.7, 0.0],
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
... ]
>>> z = mg.at_node["topographic__elevation"]
>>> z_init = mg.at_node["topographic__elevation"].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 = [
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
...     [0.0, 2.1, 1.5, 1.5, 1.6, 0.0],
...     [0.0, 2.0, 1.5, 1.5, 1.5, 0.0],
...     [0.0, 2.2, 1.5, 1.5, 1.7, 0.0],
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
... ]
>>> np.allclose(z.reshape(mg.shape), z_out)
True
>>> 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 = [
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
...     [0.0, 4.0, 8.0, 12.0, 4.0, 4.0],
...     [0.0, 4.0, 16.0, 36.0, 40.0, 40.0],
...     [0.0, 4.0, 8.0, 4.0, 4.0, 4.0],
...     [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
... ]
>>> np.allclose(mg.at_node["drainage_area"].reshape(mg.shape), 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.0  # 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
>>> lmb.lake_depths  # z was both surface and "fill_surface"
Traceback (most recent call last):
...
ValueError: surface and fill_surface must be different fields or arrays
to enable the property lake_depths!

ValueError was raised because 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.0, 4.0]))
True
>>> lmb.run_one_step()
Traceback (most recent call last):
...
ValueError: 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 because 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 behavior 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.0, 4.0]))  # 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.0)
>>> z_hex = hmg.add_zeros("topographic__elevation", at="node")
>>> z_hex[:] = hmg.node_x
>>> z_hex[11] = -3.0
>>> z_hex[12] = -1.0
>>> 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.0)
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.0)
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.0, 5.0, 3.0, 0.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).reshape(mg.shape),
...         [
...             [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
...             [0.0, 1.0, 2.0, 1.5, 4.0, 5.0, 6.0, 7.0],
...             [0.0, 1.0, 2.0, 0.5, 4.0, 5.0, 6.0, 7.0],
...             [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
...             [0.0, 1.0, 1.1, 3.0, 4.0, 5.0, 6.0, 7.0],
...             [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
...         ],
...     )
... )
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.0, m_sp=0.0, n_sp=1.0)
>>> 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"].reshape(mg.shape),
...     [
...         [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
...         [2.0, 2.0, 1.0, 4.0, 3.0, 2.0, 1.0, 0.0],
...         [1.0, 1.0, 1.0, 13.0, 3.0, 2.0, 1.0, 0.0],
...         [2.0, 2.0, 1.0, 4.0, 3.0, 2.0, 1.0, 0.0],
...         [6.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0, 0.0],
...         [0.0, 0.0, 0.0, 0.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.0)
True
>>> np.allclose(
...     mg.at_node["drainage_area"].reshape(mg.shape),
...     [
...         [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
...         [13.0, 13.0, 12.0, 4.0, 3.0, 2.0, 1.0, 0.0],
...         [2.0, 2.0, 1.0, 7.0, 3.0, 2.0, 1.0, 0.0],
...         [2.0, 2.0, 1.0, 1.0, 3.0, 2.0, 1.0, 0.0],
...         [7.0, 7.0, 6.0, 4.0, 3.0, 2.0, 1.0, 0.0],
...         [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
...     ],
... )
True
>>> np.all(
...     np.equal(
...         np.round(z_water, 2).reshape(mg.shape),
...         [
...             [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
...             [0.0, 1.0, 2.0, 2.0, 4.0, 5.0, 6.0, 7.0],
...             [0.0, 1.0, 2.0, 2.0, 4.0, 5.0, 6.0, 7.0],
...             [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
...             [0.0, 1.0, 1.1, 3.0, 4.0, 5.0, 6.0, 7.0],
...             [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
...         ],
...     )
... )
True

>>> sp.run_one_step(0.05)  # note m=0 to illustrate effect of slopes
>>> np.all(
...     np.equal(
...         np.round(z_water, 2).reshape(mg.shape),
...         [
...             [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
...             [0.0, 0.95, 1.95, 2.0, 3.9, 4.95, 5.95, 7.0],
...             [0.0, 0.95, 1.95, 2.0, 3.9, 4.95, 5.95, 7.0],
...             [0.0, 0.95, 1.95, 2.93, 3.93, 4.95, 5.95, 7.0],
...             [0.0, 0.95, 1.09, 2.91, 3.95, 4.95, 5.95, 7.0],
...             [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
...         ],
...     )
... )
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).reshape(mg.shape),
...         [
...             [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
...             [0.0, 0.95, 1.95, 1.5, 3.9, 4.95, 5.95, 7.0],
...             [0.0, 0.95, 1.95, 0.5, 3.9, 4.95, 5.95, 7.0],
...             [0.0, 0.95, 1.95, 2.93, 3.93, 4.95, 5.95, 7.0],
...             [0.0, 0.95, 1.09, 2.91, 3.95, 4.95, 5.95, 7.0],
...             [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
...         ],
...     )
... )
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).reshape(mg.shape),
...         [
...             [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
...             [0.0, 0.95, 1.95, 1.95, 3.9, 4.95, 5.95, 7.0],
...             [0.0, 0.95, 1.95, 1.95, 3.9, 4.95, 5.95, 7.0],
...             [0.0, 0.95, 1.95, 2.93, 3.93, 4.95, 5.95, 7.0],
...             [0.0, 0.95, 1.09, 2.91, 3.95, 4.95, 5.95, 7.0],
...             [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
...         ],
...     )
... )
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.

update()[source]#

Alias for running one step.

property 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, 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()
>>> lmb.was_there_overfill
Traceback (most recent call last):
...
ValueError: was_there_overfill is only defined if filling to an inclined surface!

ValueError was raised because 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,
... )
>>> lmb.run_one_step()
Traceback (most recent call last):
...
ValueError: 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 because 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, 0.2, 0.1, 1.0, 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, 0.2, 0.1, 1.0, 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.