landlab

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

Find depressions on a topographic surface.

Code author: gtucker, DEJH (Flow routing)

class DepressionFinderAndRouter(grid, routing='D8')[source]

Bases: landlab.core.model_component.Component

Find depressions on a topographic surface.

This component identifies depressions in a topographic surface, finds an outlet for each depression. If directed to do so (default True), and the component is able to find existing routing fields output from the ‘route_flow_dn’ component, it will then modify the drainage directions and accumulations already stored in the grid to route flow across these depressions.

Note that in general properties of this class named “depression” identify each individual pit in the topography, including those that will merge once the fill is performed. Those named “lake” return the unique lakes created by the fill, and are probably the properties most users will want.

Note also that the structure of drainage within the lakes is not guaranteed, and in particular, may not be symmetrical even if your boundary conditions are. However, the outputs from the lake will all still be correct.

Note the routing part of this component is not yet compatible with irregular grids.

The prinary method of this class is map_depressions(pits=’flow__sink_flag’, reroute_flow=True).

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('node', 'topographic__elevation', mg.node_x.copy())
>>> 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.

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 ('D8' or 'D4' (optional)) – 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.

__init__(grid, routing='D8')[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 ('D8' or 'D4' (optional)) – 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.

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] = CLOSED_BOUNDARY
>>> z = rg.add_zeros('node', 'topographic__elevation')
>>> z[:] = rg.x_of_node + 0.01 * rg.y_of_node
>>> lake_nodes = np.array([10, 16, 17, 18, 24, 32, 33, 38, 40])
>>> z[lake_nodes] *= 0.1
>>> fr = FlowAccumulator(rg, flow_director='D4')
>>> fr.run_one_step()
>>> rg.at_node['flow__receiver_node']
array([ 0,  1,  2,  3,  4,  5,  6,  7,  7, 16, 10, 10, 11, 13, 14, 14, 16,
       16, 17, 18, 20, 21, 21, 16, 17, 24, 33, 27, 28, 28, 29, 24, 32, 32,
       34, 35, 35, 38, 38, 38, 33, 41, 42, 43, 44, 45, 46, 47, 48])
>>> df = DepressionFinderAndRouter(rg, routing='D4')
>>> df.map_depressions()
>>> rg.at_node['flow__receiver_node']
array([ 0,  1,  2,  3,  4,  5,  6,  7,  7, 16, 17, 10, 11, 13, 14, 14, 15,
       16, 17, 18, 20, 21, 21, 16, 17, 24, 33, 27, 28, 28, 29, 38, 31, 32,
       34, 35, 35, 36, 37, 38, 33, 41, 42, 43, 44, 45, 46, 47, 48])
>>> fr = FlowAccumulator(rg, flow_director='D8')
>>> fr.run_one_step()
>>> rg.at_node['flow__receiver_node']
array([ 0,  1,  2,  3,  4,  5,  6,  7,  7, 16, 16, 10, 18, 13, 14, 14, 16,
       16, 17, 18, 20, 21, 21, 16, 16, 24, 33, 27, 28, 28, 24, 24, 24, 32,
       34, 35, 35, 38, 38, 38, 32, 41, 42, 43, 44, 45, 46, 47, 48])
>>> df = DepressionFinderAndRouter(rg, routing='D8')
>>> df.map_depressions()
>>> rg.at_node['flow__receiver_node']
array([ 0,  1,  2,  3,  4,  5,  6,  7,  7, 16, 16, 10, 18, 13, 14, 14,  8,
       16, 17, 18, 20, 21, 21, 16, 16, 24, 33, 27, 28, 28, 24, 24, 24, 32,
       34, 35, 35, 38, 32, 38, 32, 41, 42, 43, 44, 45, 46, 47, 48])
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.

find_lowest_node_on_lake_perimeter(nodes_this_depression)[source]

Locate the lowest node on the margin of the “lake”.

Parameters

nodes_this_depression (array_like of int) – Nodes that form a pit.

Returns

The lowest node on the perimeter of a depression.

Return type

int

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 LOCAL_BAD_INDEX_VALUE.

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(pits='flow__sink_flag', reroute_flow=True)[source]

Map depressions/lakes in a topographic surface.

Parameters
  • 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 ‘route_flow_dn’

  • reroute_flow (bool, optional) – If True (default), and the component detects the output fields in the grid produced by the route_flow_dn component, this component will modify the existing flow fields to route the flow across the lake surface(s). Ensure you call this method after you have already routed flow in each loop of your model.

Examples

Test #1: 5x5 raster grid with a diagonal lake.

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components.flow_routing import (
...     DepressionFinderAndRouter)
>>> rg = RasterModelGrid((5, 5))
>>> z = rg.add_zeros('node', 'topographic__elevation')
>>> z[:] = np.array([100., 100.,  95., 100., 100.,
...                  100., 101.,  92.,   1., 100.,
...                  100., 101.,   2., 101., 100.,
...                  100.,   3., 101., 101., 100.,
...                   90.,  95., 100., 100., 100.])
>>> df = DepressionFinderAndRouter(rg)
>>> df.map_depressions(pits=None, reroute_flow=False)
>>> 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.

updated_boundary_conditions()[source]

Call this if boundary conditions on the grid are updated after the component is instantiated.

PotentialityFlowRouter: Find flow directions and accumulation using potential-field theory

class PotentialityFlowRouter(self, grid, method='D8', flow_equation='default', Chezys_C=30.0, Mannings_n=0.03, **kwds)[source]

Bases: landlab.core.model_component.Component

Multidirectional flow routing using a novel method.

This class implements Voller, Hobley, and Paola’s experimental matrix solutions for flow routing. The method works by solving for a potential field at all nodes on the grid, which enforces both mass conservation and flow downhill along topographic gradients. It is order n and highly efficient, but does not return any information about flow connectivity.

Options are permitted to allow “abstract” routing (flow enforced downslope, but no particular assumptions are made about the governing equations), or routing according to the Chezy or Manning equations. This routine assumes that water is distributed evenly over the surface of the cell in deriving the depth, and does not assume channelization. You will need to back- calculate channel depths for yourself using known widths at each node if that is what you want.

It is VITAL you initialize this component AFTER setting boundary conditions.

Note that this component offers the property discharges_at_links. This returns the discharges at all links. If method==’D8’, this list will include diagonal links after the orthogonal links, which is why this information is not returned as a field.

Discharges at nodes are recorded as the outgoing total discharge (i.e., including any contribution from ‘water__unit_flux_in’).

The primary method of this class is run_one_step.

Notes

This is a “research grade” component, and is subject to dramatic change with little warning. No guarantees are made regarding its accuracy or utility. It is not recommended for user use yet!

Examples

>>> from landlab import HexModelGrid
>>> import numpy as np
>>> mg = HexModelGrid(4, 6, dx=2., node_layout='rect', orientation='vertical')
>>> z = mg.add_zeros('node', 'topographic__elevation')
>>> Q_in = mg.add_ones('node', 'water__unit_flux_in')
>>> z += mg.node_y.copy()
>>> potfr = PotentialityFlowRouter(mg)
>>> potfr.run_one_step()
>>> Q_at_core_nodes = np.array(
...     [ 13.57233404,  13.93522481,  11.52216193,  11.29307277,
...        8.80884751,   8.86380667,   6.47446459,   6.82161521])
>>> np.allclose(mg.at_node['surface_water__discharge'][mg.core_nodes],
...             Q_at_core_nodes)
True
Parameters
  • grid (ModelGrid) – A grid.

  • method ({'D8', 'D4'}, optional) – Routing method (‘D8’ is the default). This keyword has no effect for a Voronoi-based grid.

  • flow_equation ({'default', 'Manning', 'Chezy'}, optional) – If Manning or Chezy, flow is routed according to the Manning or Chezy equation; discharge is allocated to multiple downslope nodes proportional to the square root of discharge; and a water__depth field is returned. If default, flow is allocated to multiple nodes linearly with slope; and the water__depth field is not calculated.

  • Chezys_C (float (optional)) – Required if flow_equation == ‘Chezy’.

  • Mannings_n (float (optional)) – Required if flow_equation == ‘Manning’.

Return the discharges at links.

Note that if diagonal routing, this will return number_of_d8. Otherwise, it will be number_of_links.

route_flow(**kwds)[source]
run_one_step(**kwds)[source]

Route surface-water flow over a landscape.

Both convergent and divergent flow can occur.