DepressionFinderAndRouter: Handle depressions in terrain by calculating extent and drainage of “lakes”#

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

Bases: 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 may not yet be fully 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)
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)
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:

grid (RasterModelGrid) – A landlab RasterModelGrid.
routing (str) – If grid is a raster type, controls whether lake connectivity can occur on diagonals (‘D8’, default), or only orthogonally (‘D4’). Has no effect if grid is not a raster.
pits (array or str or None, optional) – If a field name, the boolean field containing True where pits. If an array, either a boolean array of nodes of the pits, or an array of pit node IDs. It does not matter whether or not open boundary nodes are flagged as pits; they are never treated as such. Default is ‘flow__sink_flag’, the pit field output from the FlowDirectors.
reroute_flow (bool, optional) – If True (default), and the component detects the output fields in the grid produced by the FlowAccumulator component, this component will modify the existing flow fields to route the flow across the lake surface(s).

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

grid (RasterModelGrid) – A landlab RasterModelGrid.
routing (str) – If grid is a raster type, controls whether lake connectivity can occur on diagonals (‘D8’, default), or only orthogonally (‘D4’). Has no effect if grid is not a raster.
pits (array or str or None, optional) – If a field name, the boolean field containing True where pits. If an array, either a boolean array of nodes of the pits, or an array of pit node IDs. It does not matter whether or not open boundary nodes are flagged as pits; they are never treated as such. Default is ‘flow__sink_flag’, the pit field output from the FlowDirectors.
reroute_flow (bool, optional) – If True (default), and the component detects the output fields in the grid produced by the FlowAccumulator component, this component will modify the existing flow fields to route the flow across the lake surface(s).

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"].reshape(rg.shape)
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"].reshape(rg.shape)
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"].reshape(rg.shape)
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"].reshape(rg.shape)
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]])

property depression_depth#: At node array of depression depths.

property depression_outlet_map#: At node array indicating the node-id of the depression outlet.

display_depression_map()[source]#: Print a simple character-based map of depressions/lakes.

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.

property flood_status#: Map of flood status (_PIT, _CURRENT_LAKE, _UNFLOODED, or _FLOODED).

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

the_node (int) – The node to test.
nodes_this_depression (array_like of int) – Nodes that form a pit.

Returns:

True if the node is a valid outlet. Otherwise, False.

Return type:

boolean

property lake_areas#

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

The order is the same as that returned by lake_codes.

property lake_at_node#: Return a boolean array, True if the node is flooded, False otherwise.

property 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”.

property 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

property lake_outlets#: Returns the unique outlets for each lake, in same order as the return from lake_codes.

property 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))
>>> rg.at_node["topographic__elevation"] = [
...     [100.0, 100.0, 95.0, 100.0, 100.0],
...     [100.0, 101.0, 92.0, 1.0, 100.0],
...     [100.0, 101.0, 2.0, 101.0, 100.0],
...     [100.0, 3.0, 101.0, 101.0, 100.0],
...     [90.0, 95.0, 100.0, 100.0, 100.0],
... ]

>>> df = DepressionFinderAndRouter(rg, reroute_flow=False)
>>> df.map_depressions()
>>> df.display_depression_map()
. . . . .
. . . ~ .
. . ~ . .
. ~ . . .
o . . . .

node_can_drain(the_node)[source]#

Check if a node has drainage away from the current lake/depression.

Parameters:

the_node (int) – The node to test.
nodes_this_depression (array_like of int) – Nodes that form a pit.

Returns:

True if the node can drain. Otherwise, False.

Return type:

boolean

property number_of_lakes#: Return the number of individual lakes.

property number_of_pits#: The number of pits on the grid.

property pit_node_ids#: Node IDs of grid nodes identified as pits.

property receivers#: At node array indicating which node receives flow.

update()[source]#: Alias for map_depressions.

updated_boundary_conditions()[source]#: Call this if boundary conditions on the grid are updated after the component is instantiated.