landlab.components.sink_fill.sink_fill_barnes¶

fill_sinks_barnes.py.

Fill sinks in a landscape to the brim, following the Barnes et al. (2014) algorithms.

class SinkFillerBarnes[source]¶

Bases: LakeMapperBarnes

Uses the Barnes et al (2014) algorithms to replace pits in a topography with flats, or optionally with very shallow gradient surfaces to allow continued draining.

This component is NOT intended for use iteratively as a model runs; rather, it is to fill in an initial topography. If you want to repeatedly fill pits as a landscape develops, you are after the LakeMapperBarnes component. If you want flow paths on your filled landscape, manually run a FlowDirector and FlowAccumulator for yourself.

The locations and depths etc. of the fills will be tracked, and properties are provided to access this information.

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

Initialise the component.

Parameters:

grid (ModelGrid) – A grid.
surface (field name at node or array of length node) – The surface to fill.
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 (at machine precision) to give steepest descent flow paths to the outlet, once they are calculated.
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.

__init__(grid, surface='topographic__elevation', method='D8', fill_flat=False, ignore_overfill=False)[source]¶

Initialise the component.

Parameters:

grid (ModelGrid) – A grid.
surface (field name at node or array of length node) – The surface to fill.
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 (at machine precision) to give steepest descent flow paths to the outlet, once they are calculated.
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.

static __new__(cls, *args, **kwds)¶

cite_as = '\n @article{BARNES2014117,\n title = {Priority-flood: An optimal depression-filling and\n watershed-labeling algorithm for digital elevation models},\n journal = "Computers & Geosciences",\n volume = "62",\n pages = "117 - 127",\n year = "2014",\n issn = "0098-3004",\n doi = "https://doi.org/10.1016/j.cageo.2013.04.024",\n url = "http://www.sciencedirect.com/science/article/pii/S0098300413001337",\n author = "Richard Barnes and Clarence Lehman and David Mulla",\n keywords = "Pit filling, Terrain analysis, Hydrology, Drainage network, Modeling, GIS"\n }'¶

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 = (('sediment_fill__depth', 'Depth of sediment added at eachnode'), ('topographic__elevation', 'Land surface topographic elevation'))¶

property fill_areas¶

A nlakes-long array of the area of each fill.

The order is the same as that of the keys in fill_dict, and of fill_outlets.

property fill_at_node¶: Return a boolean array, True if the node is filled, False otherwise.

property fill_depths¶: Return the change in surface elevation at each node this step.

property fill_dict¶

Return a dictionary where the keys are the outlet nodes of each filled area, and the values are deques of nodes within each.

Items are not returned in ID order.

property fill_map¶

Return an array of ints, where each filled node is labelled with its outlet node ID.

Nodes not in a filled area are labelled with BAD_INDEX_VALUE (default -1).

property fill_outlets¶: Returns the outlet for each filled area, not necessarily in ID order.

property fill_volumes¶

A nlakes-long array of the volume of each fill.

The order is the same as that of the keys in fill_dict, and of fill_outlets.

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.

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 = ('topographic__elevation',)¶

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

name = 'SinkFillerBarnes'¶

property number_of_fills¶: Return the number of individual filled areas.

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

optional_var_names = ()¶

output_var_names = ('sediment_fill__depth', 'topographic__elevation')¶

run_one_step()[source]¶

Fills the surface to remove all pits.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import SinkFillerBarnes, 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 = z.copy()
>>> sfb = SinkFillerBarnes(mg, method="Steepest")  # , surface=z

TODO: once return_array_at_node is fixed, this example should also take surface… GIVE IT surface=z !!

>>> sfb.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
>>> fill_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(sfb.fill_map.reshape(mg.shape) == fill_map)
True
>>> np.all(sfb.fill_at_node == (sfb.fill_map > -1))
True
>>> sfb.was_there_overfill  # everything fine with slope adding
False

>>> fr = FlowAccumulator(mg, flow_director="D4")  # routing now works
>>> fr.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, 1.0, 2.0, 3.0, 1.0, 1.0],
...     [0.0, 1.0, 4.0, 9.0, 10.0, 10.0],
...     [0.0, 1.0, 2.0, 1.0, 1.0, 1.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 and a flat fill:

>>> from collections import deque
>>> 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()
>>> sfb = SinkFillerBarnes(mg, method="Steepest", fill_flat=True)
>>> sfb.run_one_step()
>>> sfb.fill_dict == {8: deque([7]), 16: deque([15, 9, 14, 22])}
True
>>> sfb.number_of_fills
2
>>> sfb.fill_outlets == [16, 8]
True
>>> np.allclose(sfb.fill_areas, np.array([4.0, 1.0]))  # same order
True

Unlike the LakeMapperBarnes equivalents, fill_depths and fill_volumes are always available through this component:

>>> fill_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.allclose(sfb.fill_depths.reshape(mg.shape), fill_depths)
True
>>> np.allclose(sfb.fill_volumes, np.array([1.7, 1.0]))
True

Note that with a flat fill, we can’t drain the surface. The surface is completely flat, so each and every core node within the fill becomes a sink.

>>> fr.run_one_step()
>>> where_is_filled = np.where(sfb.fill_map > -1, 1, 0)
>>> np.all(
...     mg.at_node["flow__sink_flag"][mg.core_nodes]
...     == where_is_filled[mg.core_nodes]
... )
True

(Note that the fill_map does not think that the perimeter nodes are sinks, since they haven’t changed elevation. In contrast, the FlowAccumulator does think they are, because these nodes are where flow terminates.)

property shape¶: Return the grid shape attached to the component, if defined.

unit_agnostic = True¶

units = (('sediment_fill__depth', 'm'), ('topographic__elevation', 'm'))¶

update()¶: Alias for running one step.

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 = (('sediment_fill__depth', 'node'), ('topographic__elevation', '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

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.