landlab

FlowRouter: Calculate flow direction and accumulation from topography

Calculate single-path (steepest direction) flow directions.

Given a ModelGrid, calculates single-path (steepest direction) flow directions, drainage area, and (optionally) discharge.

The “dn” in the name means that this is a generalization of the D8 algorithm, for a grid in which a node has N neighbors (N might happen to be 8, or not).

class FlowRouter(*args, **kwds)[source]

Bases: landlab.core.model_component.Component

Single-path (steepest direction) flow routing.

This class implements single-path (steepest direction) flow routing, and calculates flow directions, drainage area, and (optionally) discharge.

Note that because this router is generalizd across both regular and irregular grids, perimeter nodes can NEVER contribute to the accumulating flux, even if the gradients from them point inwards to the main body of the grid. This is because under Landlab definitions, perimeter nodes lack cells, so cannot accumulate any discharge.

The primary method of this class is run_one_step().

Construction:

FlowRouter(grid, method='D8', runoff_rate=None)
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.

runoff_rate : float, optional (m/time)

If provided, sets the (spatially constant) runoff rate. If a spatially variable runoff rate is desired, use the input field ‘water__unit_flux_in’. If both the field and argument are present at the time of initialization, runoff_rate will overwrite the field. If neither are set, defaults to spatially constant unit input.

node_drainage_area
node_order_upstream
node_receiving_flow
node_steepest_slope
node_water_discharge
route_flow(**kwds)[source]

Route surface-water flow over a landscape.

Routes surface-water flow by (1) assigning to each node a single drainage direction, and then (2) adding up the number of nodes that contribute flow to each node on the grid (including the node itself).

Stores as ModelGrid fields:

  • Node array of receivers (nodes that receive flow), or ITS OWN ID if there is no receiver: ‘flow__receiver_node’
  • Node array of drainage areas: ‘drainage_area’
  • Node array of discharges: ‘surface_water__discharge’
  • Node array of steepest downhill slopes: ‘topographic__steepest_slope’
  • Node array containing downstream-to-upstream ordered list of node IDs: ‘flow__upstream_node_order’
  • Node array containing ID of link that leads from each node to its receiver, or BAD_INDEX_VALUE if no link: ‘flow__link_to_receiver_node’
  • Boolean node array of all local lows: ‘flow__sink_flag’
Returns:

ModelGrid

The modified grid object

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components.flow_routing import FlowRouter
>>> mg = RasterModelGrid((5, 4), spacing=(1, 1))
>>> elev = np.array([0.,  0.,  0., 0.,
...                  0., 21., 10., 0.,
...                  0., 31., 20., 0.,
...                  0., 32., 30., 0.,
...                  0.,  0.,  0., 0.])
>>> _ = mg.add_field('node','topographic__elevation', elev)
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fr = FlowRouter(mg)
>>> mg = fr.route_flow()
>>> mg.at_node['flow__receiver_node'] 
array([  0,  1,  2,  3,
         4,  1,  2,  7,
         8,  6,  6, 11,
        12, 10, 10, 15,
        16, 17, 18, 19])
>>> mg.at_node['drainage_area'] 
array([ 0.,  1.,  5.,  0.,
        0.,  1.,  5.,  0.,
        0.,  1.,  3.,  0.,
        0.,  1.,  1.,  0.,
        0.,  0.,  0.,  0.])

Now let’s change the cell area (100.) and the runoff rates:

>>> mg = RasterModelGrid((5, 4), spacing=(10., 10))

Put the data back into the new grid.

>>> _ = mg.add_field('node','topographic__elevation', elev)
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fr = FlowRouter(mg)
>>> runoff_rate = np.arange(mg.number_of_nodes)
>>> _ = mg.add_field('node', 'water__unit_flux_in', runoff_rate,
...                  noclobber=False)
>>> mg = fr.route_flow()
>>> mg.at_node['surface_water__discharge'] 
array([    0.,   500.,  5200.,     0.,
           0.,   500.,  5200.,     0.,
           0.,   900.,  3700.,     0.,
           0.,  1300.,  1400.,     0.,
           0.,     0.,     0.,     0.])
run_one_step(**kwds)[source]

Route surface-water flow over a landscape.

Routes surface-water flow by (1) assigning to each node a single drainage direction, and then (2) adding up the number of nodes that contribute flow to each node on the grid (including the node itself).

This is the fully standardized run method for this component. It differs from route_flow() in that it has a standardized name, and does not return anything.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components.flow_routing import FlowRouter
>>> mg = RasterModelGrid((5, 4), spacing=(1, 1))
>>> elev = np.array([0.,  0.,  0., 0.,
...                  0., 21., 10., 0.,
...                  0., 31., 20., 0.,
...                  0., 32., 30., 0.,
...                  0.,  0.,  0., 0.])
>>> _ = mg.add_field('node','topographic__elevation', elev)
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fr = FlowRouter(mg)
>>> fr.run_one_step()
>>> mg.at_node['flow__receiver_node'] 
array([  0,  1,  2,  3,
         4,  1,  2,  7,
         8,  6,  6, 11,
        12, 10, 10, 15,
        16, 17, 18, 19])
>>> mg.at_node['drainage_area'] 
array([ 0.,  1.,  5.,  0.,
        0.,  1.,  5.,  0.,
        0.,  1.,  3.,  0.,
        0.,  1.,  1.,  0.,
        0.,  0.,  0.,  0.])

Now let’s change the cell area (100.) and the runoff rates:

>>> mg = RasterModelGrid((5, 4), spacing=(10., 10))

Put the data back into the new grid.

>>> _ = mg.add_field('node','topographic__elevation', elev)
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fr = FlowRouter(mg)
>>> runoff_rate = np.arange(mg.number_of_nodes)
>>> _ = mg.add_field('node', 'water__unit_flux_in', runoff_rate,
...                  noclobber=False)
>>> fr.run_one_step()
>>> mg.at_node['surface_water__discharge'] 
array([    0.,   500.,  5200.,     0.,
           0.,   500.,  5200.,     0.,
           0.,   900.,  3700.,     0.,
           0.,  1300.,  1400.,     0.,
           0.,     0.,     0.,     0.])
updated_boundary_conditions()[source]

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

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', **kwds)[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).

Construction:

DepressionFinderAndRouter(grid, grid, routing='D8')
Parameters:

grid : RasterModelGrid

A landlab RasterModelGrid.

routing : {‘D8’, ‘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.

Examples

Route flow across a depression in a sloped surface.

>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowRouter, DepressionFinderAndRouter
>>> mg = RasterModelGrid((7, 7), 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 = FlowRouter(mg)
>>> 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 FlowRouter will have now been modified to route flow over the depressions in the surface. The topogrphy itself is not modified.

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 FlowRouter
>>> 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 = FlowRouter(rg, method='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 = FlowRouter(rg, method='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:

int

The lowest node on the perimeter of a depression.

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:

boolean

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

lake_areas

A nlakes-long array of the area of each lake. The order is the same as that returned by lake_codes.

lake_at_node

Return a boolean array, True if the node is flooded, False otherwise.

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

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.

lake_outlets

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

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:

boolean

True if the node can drain. Otherwise, False.

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.

main()[source]

temporary: test.

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

class PotentialityFlowRouter(*args, **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().

Construction:

PotentialityFlowRouter(grid, method='D8', flow_equation='default',
             Chezys_C=30., Mannings_n=0.03)
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’.

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., shape='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(
...     [ 17.02012846,  16.88791903,  13.65746194,  14.85578934,
...       11.41908145,  11.43630865,   8.95902559,  10.04348075,
...        6.28696459,   6.44316089,   4.62478522,   5.29145188])
>>> np.allclose(mg.at_node['surface_water__discharge'][mg.core_nodes],
...             Q_at_core_nodes)
True

Return the discharges at links.

Note that if diagonal routing, this will return number_of_d8_links. 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.