# The Landlab FlowDirectors: Components for Flow Direction¶

## FlowDirectorSteepest¶

flow_director_steepest.py: provides the component FlowDirectorSteepest.

This components finds the steepest single-path steepest descent flow directions. It is equivalent to D4 method in the special case of a raster grid in that it does not consider diagonal links between nodes. For that capability, use FlowDirectorD8.

class FlowDirectorSteepest(grid, surface='topographic__elevation')[source]

Bases: landlab.components.flow_director.flow_director_to_one._FlowDirectorToOne

Single-path (steepest direction) flow direction without diagonals.

This components finds the steepest single-path steepest descent flow directions. It is equivalent to D4 method in the special case of a raster grid in that it does not consider diagonal links between nodes. For that capability, use FlowDirectorD8.

Stores as ModelGrid fields:

• Node array of steepest downhill slopes: ‘topographic__steepest_slope’
• Boolean node array of all local lows: ‘flow__sink_flag’

The primary method of this class is run_one_step.

Examples

This method works for both raster and irregular grids. First we will look at a raster example, and then an irregular example.

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorSteepest
>>> mg = RasterModelGrid((3,3), xy_spacing=(1, 1))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
...     "topographic__elevation",
...     mg.node_x + mg.node_y,
...     at="node",
... )
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fd.surface_values
array([ 0.,  1.,  2.,  1.,  2.,  3.,  2.,  3.,  4.])
>>> fd.run_one_step()
array([0, 1, 2, 3, 1, 5, 6, 7, 8])
>>> mg.at_node['topographic__steepest_slope']
array([ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.])
array([-1, -1, -1, -1,  3, -1, -1, -1, -1])
>>> mg.at_node['flow__sink_flag'].astype(int)
array([1, 1, 1, 1, 0, 1, 1, 1, 1])
>>> mg_2 = RasterModelGrid((5, 4), xy_spacing=(1, 1))
>>> topographic__elevation = np.array([0.,  0.,  0., 0.,
...                                    0., 21., 10., 0.,
...                                    0., 31., 20., 0.,
...                                    0., 32., 30., 0.,
...                                    0.,  0.,  0., 0.])
...     "topographic__elevation",
...     topographic__elevation,
...     at="node",
... )
>>> mg_2.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fd_2 = FlowDirectorSteepest(mg_2)
>>> fd_2.run_one_step()
array([ 0,  1,  2,  3,
4,  1,  2,  7,
8, 10,  6, 11,
12, 14, 10, 15,
16, 17, 18, 19])


And the at-link field 'flow__link_direction' indicates if the flow along the link is with or against the direction indicated by 'link_dirs_at_node' (from tail node to head node).

>>> mg_2.at_link['flow__link_direction']
array([ 0,  0,  0,  0, -1, -1,  0,  0,  0,  0,  0,  0, -1,  0,  0,  1,  0,
0,  0, -1,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0], dtype=int8)


This indicates that flow on links 4, 5, 12, and 19 goes against the topologic ordering – that is that flow goes from head node to tail node – and that flow goes with the topologic ordering on links 15 and 22. All other links have no flow on them.

The FlowDirectorSteepest attribute flow_link_direction_at_node indicates the link flow direction (with or against topology directions) for all links at node. The ordering of links at node mirrors the grid attribute links_at_node.

>>> fd_2.flow_link_direction_at_node()
array([[ 0,  0,  0,  0],
[ 0, -1,  0,  0],
[ 0, -1,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0, -1],
[ 0, -1,  0, -1],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 1,  0,  0,  0],
[ 0, -1,  1, -1],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 1,  0,  0,  0],
[ 0,  0,  1, -1],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0]], dtype=int8)


For example, this indicates that node 10 has flow going along three links that are attached to it. The link to the East has no flow, the link to the North has flow going against the topologic direction, the link to the West has flow going with the topologic direction, and the link to the South has flow going against the topologic direction.

In many use cases, one might want to know which links are bringing flow into or out of the node. The flow director attribute flow_link_incoming_at_node provides this information. Here -1 means that flow is outgoing from the node and 1 means it is incoming.

>>> fd_2.flow_link_incoming_at_node()
array([[ 0,  0,  0,  0],
[ 0,  1,  0,  0],
[ 0,  1,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0, -1],
[ 0,  1,  0, -1],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[-1,  0,  0,  0],
[ 0,  1,  1, -1],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[-1,  0,  0,  0],
[ 0,  0,  1, -1],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0]], dtype=int8)


So if one wanted to identify the source nodes at node, you would do the following:

>>> np.where(fd_2.flow_link_incoming_at_node() == 1, mg_2.adjacent_nodes_at_node, -1)
array([[-1, -1, -1, -1],
[-1,  5, -1, -1],
[-1,  6, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, 10, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, 14,  9, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, -1, 13, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1],
[-1, -1, -1, -1]])


The flow directors also have the ability to return the flow receiver nodes

>>> receiver = fd.direct_flow()
array([0, 1, 2,
3, 1, 5,
6, 7, 8])


For the second example we will use a Hexagonal Model Grid, a special type of Voroni Grid that has regularly spaced hexagonal cells.

>>> from landlab import HexModelGrid
>>> mg = HexModelGrid((5, 3))
...     "topographic__elevation",
...     mg.node_x + np.round(mg.node_y),
...     at="node",
... )
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fd.surface_values
array([ 1. ,  2. ,  3. ,
1.5,  2.5,  3.5,  4.5,
2. ,  3. ,  4. ,  5. ,  6. ,
3.5,  4.5,  5.5,  6.5,
4. ,  5. ,  6. ])
>>> fd.run_one_step()
array([ 0,  1,  2,
3,  0,  1,  6,
7,  3,  4,  5,  11,
12,  8,  9, 15,
16, 17, 18])
>>> mg.at_node['topographic__steepest_slope']
array([ 0. ,  0. ,  0. ,
0. ,  1.5,  1.5,   0. ,
0. ,  1.5,  1.5,  1.5,  0. ,
0. ,  1.5,  1.5,  0. ,
0. ,  0. ,  0. ])
array([-1, -1, -1,
-1,  3,  5, -1,
-1, 12, 14, 16, -1,
-1, 25, 27, -1,
-1, -1, -1])
>>> mg.at_node['flow__sink_flag'].astype(int)
array([1, 1, 1,
1, 0, 0, 1,
1, 0, 0, 0, 1,
1, 0, 0, 1,
1, 1, 1])
array([ 0,  1,  2,
3,  0,  1,  6,
7,  3,  4,  5, 11,
12,  8,  9, 15,
16, 17, 18])


References

Required Software Citation(s) Specific to this Component

None Listed

None Listed

Parameters: grid (ModelGrid) – A grid. surface (field name at node or array of length node, optional) – The surface to direct flow across, default is field at node: topographic__elevation,.
__init__(grid, surface='topographic__elevation')[source]
Parameters: grid (ModelGrid) – A grid. surface (field name at node or array of length node, optional) – The surface to direct flow across, default is field at node: topographic__elevation,.
direct_flow()[source]

Find flow directions, save to the model grid, and return receivers.

direct_flow() checks for updated boundary conditions, calculates slopes on links, finds baselevel nodes based on the status at node, calculates flow directions, saves results to the grid, and returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_node’]

An alternative to direct_flow() is run_one_step() which does the same things but also returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_node’]

At-link array of the downstream node based on flow direction.

BAD_INDEX_VALUE is given if no downstream node is defined.

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorSteepest
>>> mg = RasterModelGrid((3,3))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
...     "topographic__elevation",
...     mg.node_x + mg.node_y,
...     at="node",
... )
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fd.run_one_step()
array([-1, -1, -1,  1, -1, -1, -1, -1, -1, -1, -1, -1])


Return array of flow link direction.

This property indicates the relationship between the flow direction (determined based on the elevation of nodes) and the topologic link direction (in which the head and tail nodes are defined based on relative position in x-y space).

Recall that the standard landlab link direction goes from the tail node to the head node.

A value of zero indicates that the link does not exist or is not active.

A value of -1 indicates that water flow based on flow__link_to_receiver_node goes from head node to tail node, while a value of 1 indicates that water flow goes from tail node to head node.

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorSteepest
>>> mg = RasterModelGrid((3,3))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
...     "topographic__elevation",
...     mg.node_x + mg.node_y,
...     at="node",
... )
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fd.run_one_step()
array([ 0,  0,  0, -1,  0,  0,  0,  0,  0,  0,  0,  0], dtype=int8)


Return array of flow link direction at node.

This property mirrors links_at_node and indicates the relationship between the flow direction (determined based on the elevation of nodes) and the topologic link direction (in which the head and tail nodes are defined based on relative position in x-y space).

It has the shape (number of nodes, maximum number of links at node).

Recall that the standard landlab link direction goes from the tail node to the head node.

A value of zero indicates that the link does not exist or is not active.

A value of -1 indicates that water flow based on flow__link_to_receiver_node goes from head node to tail node, while a value of 1 indicates that water flow goes from tail node to head node.

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorSteepest
>>> mg = RasterModelGrid((3,3))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
...     "topographic__elevation",
...     mg.node_x + mg.node_y,
...     at="node",
... )
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fd.run_one_step()
array([[ 0,  0,  0,  0],
[ 0, -1,  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,  0,  0]], dtype=int8)


This method will be updated when the DepressionFinderAndRouter is run.

First, without DepressionFinderAndRouter:

>>> from landlab.components import FlowAccumulator
>>> mg1 = RasterModelGrid((5, 5))
...     "topographic__elevation",
...     mg1.x_of_node+2 * mg1.y_of_node,
...     at="node",
... )
>>> z1[12] -= 5
>>> mg1.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fa1 = FlowAccumulator(mg1, flow_director='Steepest')
>>> fa1.run_one_step()
array([-1, -1, -1, -1, -1,
-1,  5, 15,  7, -1,
-1, 19, -1, 20, -1,
-1, 23, 24, 25, -1,
-1, -1, -1, -1, -1])
array([[ 0,  0,  0,  0],
[ 0, -1,  0,  0],
[ 0,  0,  0,  0],
[ 0, -1,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0, -1],
[ 0,  1,  0,  0],
[ 0,  0,  0, -1],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 1, -1,  0,  0],
[-1, -1,  1,  1],
[ 0, -1, -1,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0, -1],
[ 0,  0,  0, -1],
[ 0,  0,  0, -1],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0]], dtype=int8)


Next with DepressionFinderAndRouter:

>>> mg2 = RasterModelGrid((5, 5))
...     "topographic__elevation",
...     mg2.x_of_node+2 * mg2.y_of_node,
...     at="node",
... )
>>> z2[12] -= 5
>>> mg2.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fa2 = FlowAccumulator(
...       mg2,
...       flow_director='Steepest',
...       depression_finder='DepressionFinderAndRouter',
...       routing='D4'
... )
>>> fa2.run_one_step()
array([-1, -1, -1, -1, -1,
-1,  5,  6,  7, -1,
-1, 19, 15, 20, -1,
-1, 23, 24, 25, -1,
-1, -1, -1, -1, -1])
array([[ 0,  0,  0,  0],
[ 0, -1,  0,  0],
[ 0, -1,  0,  0],
[ 0, -1,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0, -1],
[ 0, -1,  0, -1],
[ 0,  0,  0, -1],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 1, -1,  0,  0],
[-1, -1,  1, -1],
[ 0, -1, -1,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0, -1],
[ 0,  0,  0, -1],
[ 0,  0,  0, -1],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0],
[ 0,  0,  0,  0]], dtype=int8)


Return array that mirrors links at node and indicates incoming flow.

This array has the shape (number of nodes, maximum number of links at node).

Incoming flow is indicated as 1 and outgoing as -1. 0 indicates that no flow moves along the link or that the link does not exist.

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorSteepest
>>> mg = RasterModelGrid((3,3))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
...     "topographic__elevation",
...     mg.node_x + mg.node_y,
...     at="node",
... )
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fd.run_one_step()
array([[ 0,  0,  0,  0],
[ 0,  1,  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,  0,  0]], dtype=int8)

run_one_step()[source]

Find flow directions and save to the model grid.

run_one_step() checks for updated boundary conditions, calculates slopes on links, finds baselevel nodes based on the status at node, calculates flow directions, and saves results to the grid.

An alternative to direct_flow() is run_one_step() which does the same things but also returns the receiver nodes not return values.

updated_boundary_conditions()[source]

Method to update FlowDirectorSteepest when boundary conditions change.

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

At-link array of the upstream node based on flow direction.

BAD_INDEX_VALUE is given if no upstream node is defined.

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorSteepest
>>> mg = RasterModelGrid((3,3))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
...     "topographic__elevation",
...     mg.node_x + mg.node_y,
...     at="node",
... )
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fd.run_one_step()
array([-1, -1, -1,  4, -1, -1, -1, -1, -1, -1, -1, -1])


## FlowDirectorD8¶

flow_director_d8.py: provides the component FlowDirectorsD8.

This components finds the steepest single-path steepest descent flow directions and considers diagonal links between nodes on a raster grid. It is not implemented for irregular grids. For a method that works for irregular grids and does not consider diagonal links for rasters, use FlowDirectorSteepest instead.

class FlowDirectorD8(grid, surface='topographic__elevation')[source]

Bases: landlab.components.flow_director.flow_director_to_one._FlowDirectorToOne

Single-path (steepest direction) flow direction with diagonals on rasters.

Single-path (steepest direction) flow direction finding on raster grids by the D8 method. This method considers flow on all eight links such that flow is possible on orthogonal and on diagonal links.

The method that considers only orthogonal links (D4 method) for raster grids is FlowDirectorSteepest.

This method is not implemented for Voroni grids, use FlowDirectorSteepest instead.

Stores as ModelGrid fields:

• Node array of steepest downhill slopes: ‘topographic__steepest_slope’
• Boolean node array of all local lows: ‘flow__sink_flag’

The primary method of this class is run_one_step.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorD8
>>> mg = RasterModelGrid((3,3), xy_spacing=(1, 1))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
...     "topographic__elevation",
...     mg.node_x + mg.node_y,
...     at="node",
... )
>>> fd = FlowDirectorD8(mg, 'topographic__elevation')
>>> fd.surface_values
array([ 0.,  1.,  2.,  1.,  2.,  3.,  2.,  3.,  4.])
>>> fd.run_one_step()
array([0, 1, 2, 3, 0, 5, 6, 7, 8])
>>> mg.at_node['topographic__steepest_slope']
array([ 0.        ,  0.        ,  0.        ,  0.        ,  1.41421356,
0.        ,  0.        ,  0.        ,  0.        ])
array([-1, -1, -1, -1, 12, -1, -1, -1, -1])
>>> mg.at_node['flow__sink_flag'].astype(int)
array([1, 1, 1, 1, 0, 1, 1, 1, 1])
>>> mg_2 = RasterModelGrid((5, 4), xy_spacing=(1, 1))
>>> topographic__elevation = np.array([0.,  0.,  0., 0.,
...                                    0., 21., 10., 0.,
...                                    0., 31., 20., 0.,
...                                    0., 32., 30., 0.,
...                                    0.,  0.,  0., 0.])
...     "topographic__elevation",
...     topographic__elevation,
...     at="node",
... )
>>> mg_2.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fd_2 = FlowDirectorD8(mg_2)
>>> fd_2.run_one_step()
array([  0,  1,  2,  3,
4,  1,  2,  7,
8,  6,  6, 11,
12, 10, 10, 15,
16, 17, 18, 19])


The flow directors also have the ability to return the flow receiver nodes

>>> receiver = fd.direct_flow()
array([0, 1, 2,
3, 0, 5,
6, 7, 8])


References

Required Software Citation(s) Specific to this Component

None Listed

O’Callaghan, J., Mark, D. (1984). The extraction of drainage networks from digital elevation data. Computer Vision, Graphics, and Image Processing 28(3), 323 - 344. https://dx.doi.org/10.1016/s0734-189x(84)80011-0

Parameters: grid (ModelGrid) – A grid of type RasterModelGrid. surface (field name at node or array of length node, optional) – The surface to direct flow across, default is field at node: topographic__elevation,.
__init__(grid, surface='topographic__elevation')[source]
Parameters: grid (ModelGrid) – A grid of type RasterModelGrid. surface (field name at node or array of length node, optional) – The surface to direct flow across, default is field at node: topographic__elevation,.
direct_flow()[source]

Find flow directions, save to the model grid, and return receivers.

direct_flow() checks for updated boundary conditions, calculates slopes on links, finds baselevel nodes based on the status at node, calculates flow directions, saves results to the grid, and returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_node’]

an alternative to direct_flow() is run_one_step() which does the same things but also returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_node’]

run_one_step()[source]

Find flow directions and save to the model grid.

run_one_step() checks for updated boundary conditions, calculates slopes on links, finds baselevel nodes based on the status at node, calculates flow directions, and saves results to the grid.

an alternative to direct_flow() is direct_flow() which does the same things but also returns the receiver nodes not return values.

updated_boundary_conditions()[source]

Method to update FlowDirectorD8 when boundary conditions change.

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

## FlowDirectorMFD¶

flow_director_mfd.py: provides the component FlowDirectorMFD.

This components finds the steepest single-path steepest descent flow directions. It is equivalent to D4 method in the special case of a raster grid in that it does not consider diagonal links between nodes. For that capability, use FlowDirectorD8.

class FlowDirectorMFD(grid, surface='topographic__elevation', **kwargs)[source]

Bases: landlab.components.flow_director.flow_director_to_many._FlowDirectorToMany

Multiple-path flow direction with or without out diagonals.

Directs flow by the multiple flow direction method. Each node is assigned multiple flow directions, toward all of the N neighboring nodes that are lower than it. If none of the neighboring nodes are lower, the location is identified as a pit. Flow proportions can be calculated as proportional to slope or proportional to the square root of slope, which is the solution to a steady kinematic wave.

Specifically, it stores as ModelGrid fields:

• Node array of receivers (nodes that receive flow), or ITS OWN ID if there is no receiver: ‘flow__receiver_node’. This array is 2D, and is of dimension (number of nodes x max number of receivers).
• Node array of flow proportions: ‘flow__receiver_proportions’. This array is 2D, and is of dimension (number of nodes x max number of receivers).
• Node array of links carrying flow: ‘flow__link_to_receiver_node’. This array is 2D, and is of dimension (number of nodes x max number of receivers).
• Node array of downhill slopes corresponding to each receiver.: ‘topographic__steepest_slope’ This array is 2D, and is of dimension (number of nodes x max number of receivers). If the slope is uphill or flat, the value is assigned zero.
• Boolean node array of all local lows: ‘flow__sink_flag’

The primary method of this class is run_one_step.

Examples

This method works for both raster and irregular grids. First we will look at a raster example, and then an irregular example.

>>> import numpy as numpy
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorMFD
>>> mg = RasterModelGrid((3,3), xy_spacing=(1, 1))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
...     "topographic__elevation",
...     mg.node_x + mg.node_y,
...     at="node",
... )


The MFD flow director can be uses for raster and irregular grids. For raster grids, use of diagonal links is specified with the keyword diagonals (default is False).

>>> fd = FlowDirectorMFD(mg, 'topographic__elevation', diagonals = True)
>>> fd.surface_values
array([ 0.,  1.,  2.,  1.,  2.,  3.,  2.,  3.,  4.])
>>> fd.run_one_step()


Unlike flow directors that only direct flow to one node, FlowDirectorMFD directs flow to all downstream nodes. It stores the receiver information is a (number of nodes x maximum number or receivers) shape field at nodes.

>>> mg.at_node['flow__receiver_node']
array([[ 0, -1, -1, -1, -1, -1, -1, -1],
[ 1, -1, -1, -1, -1, -1, -1, -1],
[ 2, -1, -1, -1, -1, -1, -1, -1],
[ 3, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1,  1, -1, -1,  0, -1],
[ 5, -1, -1, -1, -1, -1, -1, -1],
[ 6, -1, -1, -1, -1, -1, -1, -1],
[ 7, -1, -1, -1, -1, -1, -1, -1],
[ 8, -1, -1, -1, -1, -1, -1, -1]])


It also stores the proportions of flow going to each receiver, the link on which the flow moves in at node arrays, and the slope of each link.

>>> mg.at_node['flow__receiver_proportions'] # doctest: +NORMALIZE_WHITESPACE
array([[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ],
[ 0.        ,  0.        ,  0.        ,  0.41421356,  0.        ,
0.        ,  0.58578644,  0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ]])
array([[-1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1,  3, -1, -1, 12, -1],
[-1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1, -1]])
>>> mg.at_node['topographic__steepest_slope'] # doctest: +NORMALIZE_WHITESPACE
array([[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ],
[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ],
[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ],
[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ],
[ 0.        ,  0.        ,  0.        ,  1.        ,  0.        ,
0.        ,  1.41421356,  0.        ],
[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ],
[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ],
[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ],
[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ]])


Finally, FlowDirectorMFD identifies sinks, or local lows.

>>> mg.at_node['flow__sink_flag'].astype(int)
array([1, 1, 1, 1, 0, 1, 1, 1, 1])


The flow directors also have the ability to return the flow receiver nodes. For this example, we will turn the diagonals off. This is the default value.

>>> mg = RasterModelGrid((3,3), xy_spacing=(1, 1))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
...     "topographic__elevation",
...     mg.node_x + mg.node_y,
...     at="node",
... )
>>> fd = FlowDirectorMFD(mg, 'topographic__elevation')
>>> fd.run_one_step()
array([[ 0, -1, -1, -1],
[ 1, -1, -1, -1],
[ 2, -1, -1, -1],
[ 3, -1, -1, -1],
[-1, -1, -1,  1],
[ 5, -1, -1, -1],
[ 6, -1, -1, -1],
[ 7, -1, -1, -1],
[ 8, -1, -1, -1]])
>>> proportions # doctest: +NORMALIZE_WHITESPACE
array([[ 1.,  0.,  0.,  0.],
[ 1.,  0.,  0.,  0.],
[ 1.,  0.,  0.,  0.],
[ 1.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  1.],
[ 1.,  0.,  0.,  0.],
[ 1.,  0.,  0.,  0.],
[ 1.,  0.,  0.,  0.],
[ 1.,  0.,  0.,  0.]])


For each donor node (represented by each row) the proportions should sum to one.

>>> proportions.sum(axis=1)
array([ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.])


For the second example we will use a Hexagonal Model Grid, a special type of Voroni Grid that has regularly spaced hexagonal cells. FlowDirectorMFD has multiple ways to partition flow based on slope. The default method is based on the slope angle. A secondary methods is to partion based on the square root of slope. This represents the solution to a steady kinematic wave.

>>> from landlab import HexModelGrid
>>> mg = HexModelGrid((5, 3))
...     "topographic__elevation",
...     mg.node_x + numpy.round(mg.node_y),
...     at="node",
... )
>>> fd = FlowDirectorMFD(
...      mg,
...      'topographic__elevation',
...      partition_method='square_root_of_slope'
...      )
>>> fd.surface_values # doctest: +NORMALIZE_WHITESPACE
array([ 1. ,  2. ,  3. ,
1.5,  2.5,  3.5,  4.5,
2. ,  3. ,  4. ,  5. ,  6. ,
3.5,  4.5,  5.5,  6.5,
4. ,  5. ,  6. ])
>>> fd.run_one_step()
array([[ 0, -1, -1, -1, -1, -1],
[ 1, -1, -1, -1, -1, -1],
[ 2, -1, -1, -1, -1, -1],
[ 3, -1, -1, -1, -1, -1],
[-1, -1, -1,  3,  0,  1],
[-1, -1, -1,  4,  1,  2],
[ 6, -1, -1, -1, -1, -1],
[ 7, -1, -1, -1, -1, -1],
[-1, -1, -1,  7,  3,  4],
[-1, -1, -1,  8,  4,  5],
[-1, -1, -1,  9,  5,  6],
[11, -1, -1, -1, -1, -1],
[12, -1, -1, -1, -1, -1],
[-1, -1, 16, 12,  8,  9],
[-1, -1, 17, 13,  9, 10],
[15, -1, -1, -1, -1, -1],
[16, -1, -1, -1, -1, -1],
[17, -1, -1, -1, -1, -1],
[18, -1, -1, -1, -1, -1]])
array([[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ],
[ 0.        ,  0.        ,  0.        ,  0.34108138,  0.41773767,
0.24118095],
[ 0.        ,  0.        ,  0.        ,  0.34108138,  0.41773767,
0.24118095],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ],
[ 0.        ,  0.        ,  0.        ,  0.34108138,  0.41773767,
0.24118095],
[ 0.        ,  0.        ,  0.        ,  0.34108138,  0.41773767,
0.24118095],
[ 0.        ,  0.        ,  0.        ,  0.34108138,  0.41773767,
0.24118095],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ],
[ 0.        ,  0.        ,  0.19431571,  0.27480391,  0.33656468,
0.19431571],
[ 0.        ,  0.        ,  0.19431571,  0.27480391,  0.33656468,
0.19431571],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ]])
array([[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1,  8,  3,  4],
[-1, -1, -1,  9,  5,  6],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, 19, 12, 13],
[-1, -1, -1, 20, 14, 15],
[-1, -1, -1, 21, 16, 17],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1],
[-1, -1, 35, 31, 25, 26],
[-1, -1, 37, 32, 27, 28],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1]])
>>> mg.at_node['topographic__steepest_slope'] # doctest: +NORMALIZE_WHITESPACE
array([[ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ],
[ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ],
[ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ],
[ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ],
[ 0. ,  0. ,  0. ,  1. ,  1.5,  0.5],
[ 0. ,  0. ,  0. ,  1. ,  1.5,  0.5],
[ 0. ,  0. ,  1. ,  0. ,  0. ,  0. ],
[ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ],
[ 0. ,  0. ,  0. ,  1. ,  1.5,  0.5],
[ 0. ,  0. ,  0. ,  1. ,  1.5,  0.5],
[ 0. ,  0. ,  0. ,  1. ,  1.5,  0.5],
[ 0. ,  1. ,  0. ,  0. ,  0. ,  0. ],
[ 0. ,  0. ,  0. ,  0.5,  0. ,  0. ],
[ 0. ,  0. ,  0.5,  1. ,  1.5,  0.5],
[ 0. ,  0. ,  0.5,  1. ,  1.5,  0.5],
[ 0. ,  1. ,  1.5,  0. ,  0. ,  0. ],
[ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ],
[ 0. ,  0. ,  0.5,  0. ,  0. ,  0. ],
[ 0. ,  0.5,  0. ,  0. ,  0. ,  0. ]])
>>> mg.at_node['flow__sink_flag'].astype(int)
array([1, 1, 1,
1, 0, 0, 1,
1, 0, 0, 0, 1,
1, 0, 0, 1,
1, 1, 1])


References

Required Software Citation(s) Specific to this Component

None Listed

Freeman, T. (1991). Calculating catchment area with divergent flow based on a regular grid. Computers and Geosciences 17(3), 413 - 422. https://dx.doi.org/10.1016/0098-3004(91)90048-i

Quinn, P., Beven, K., Chevallier, P., Planchon, O. (1991). The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models Hydrological Processes 5(1), 59-79. https://dx.doi.org/10.1002/hyp.3360050106

Parameters: grid (ModelGrid) – A grid. surface (field name at node or array of length node, optional) – The surface to direct flow across, default is field at node: topographic__elevation. partition_method (string, optional) – Method for partitioning flow. Options include ‘slope’ (default) and ‘square_root_of_slope’.
__init__(grid, surface='topographic__elevation', **kwargs)[source]
Parameters: grid (ModelGrid) – A grid. surface (field name at node or array of length node, optional) – The surface to direct flow across, default is field at node: topographic__elevation. partition_method (string, optional) – Method for partitioning flow. Options include ‘slope’ (default) and ‘square_root_of_slope’.
direct_flow()[source]

Find flow directions, save to the model grid, and return receivers.

direct_flow() checks for updated boundary conditions, calculates slopes on links, finds basself._surface_valuesel nodes based on the status at node, calculates flow directions, saves results to the grid, and returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_nodes’]

An alternative to direct_flow() is run_one_step() which does the same things but also returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_nodes’]

run_one_step()[source]

Find flow directions and save to the model grid.

run_one_step() checks for updated boundary conditions, calculates slopes on links, finds basself._surface_valuesel nodes based on the status at node, calculates flow directions, and saves results to the grid.

An alternative to run_one_step() is direct_flow() which does the same things but also returns the receiver nodes not return values.

updated_boundary_conditions()[source]

Method to update FlowDirectorMFD when boundary conditions change.

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

## FlowDirectorDinf¶

flow_director_dinf.py: provides the component FlowDirectorDINF.

Directs flow on raster grids only using the Dinfinity algorithm of Tarboton 1997.

class FlowDirectorDINF(grid, surface='topographic__elevation')[source]

Bases: landlab.components.flow_director.flow_director_to_many._FlowDirectorToMany

Flow direction on a raster grid by the D infinity method.

Directs flow by the D infinity method (Tarboton, 1997). Each node is assigned two flow directions, toward the two neighboring nodes that are on the steepest subtriangle. Partitioning of flow is done based on the aspect of the subtriangle.

Specifically, it stores as ModelGrid fields:

• Node array of receivers (nodes that receive flow), or ITS OWN ID if there is no receiver: ‘flow__receiver_node’. This array is 2D, and is of dimension (number of nodes x max number of receivers).
• Node array of flow proportions: ‘flow__receiver_proportions’. This array is 2D, and is of dimension (number of nodes x max number of receivers).
• Node array of links carrying flow: ‘flow__link_to_receiver_node’. This array is 2D, and is of dimension (number of nodes x max number of receivers).
• Node array of downhill slopes corresponding to each receiver.: ‘topographic__steepest_slope’ This array is 2D, and is of dimension (number of nodes x max number of receivers). If the slope is uphill or flat, the value is assigned zero.
• Boolean node array of all local lows: ‘flow__sink_flag’

The primary method of this class is run_one_step.

Examples

This method works for both raster and irregular grids. First we will look at a raster example, and then an irregular example.

>>> import numpy as numpy
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowDirectorDINF
>>> mg = RasterModelGrid((4,4), xy_spacing=(1, 1))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
...     "topographic__elevation",
...     mg.node_x**2 + mg.node_y**2,
...     at="node",
... )


The DINF flow director can be uses for raster grids only.

>>> fd = FlowDirectorDINF(mg, 'topographic__elevation')
>>> fd.surface_values.reshape(mg.shape)
array([[  0.,   1.,   4.,   9.],
[  1.,   2.,   5.,  10.],
[  4.,   5.,   8.,  13.],
[  9.,  10.,  13.,  18.]])
>>> fd.run_one_step()


Unlike flow directors that only direct flow to one node or to all downstream nodes, FlowDirectorDINF directs flow two nodes only. It stores the receiver information is a (number of nodes x 2) shape field at nodes.

>>> mg.at_node['flow__receiver_node']
array([[ 0, -1],
[ 1, -1],
[ 2, -1],
[ 3, -1],
[ 4, -1],
[ 0,  -1],
[ 5,  1],
[ 7, -1],
[ 8, -1],
[ 5, -1],
[ 5, -1],
[11, -1],
[12, -1],
[13, -1],
[14, -1],
[15, -1]])


It also stores the proportions of flow going to each receiver, the link on which the flow moves in at node arrays, and the slope of each link.

>>> mg.at_node['flow__receiver_proportions'] # doctest: +NORMALIZE_WHITESPACE
array([[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 0.59033447,  0.40966553],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ]])
array([[-1, -1],
[-1, -1],
[-1, -1],
[-1, -1],
[ 3, 25],
[24,  4],
[ 8, 26],
[ 9, 28],
[14, 31],
[11, 33],
[32, 12],
[16, 34],
[21, 37],
[18, 39],
[19, 38],
[20, 40]])
>>> mg.at_node['topographic__steepest_slope'] # doctest: +NORMALIZE_WHITESPACE
array([[ -1.00000000e+00,  -1.41421356e+00],
[  1.00000000e+00,  -7.12763635e+02],
[  3.00000000e+00,   1.41421356e+00],
[  5.00000000e+00,   2.82842712e+00],
[  1.00900000e+03,   7.12763635e+02],
[  1.41421356e+00,   1.00000000e+00],
[  3.00000000e+00,   2.82842712e+00],
[  1.00400000e+03,   7.10642315e+02],
[  1.00400000e+03,   7.12056529e+02],
[  3.00000000e+00,   0.00000000e+00],
[  4.24264069e+00,   3.00000000e+00],
[  1.00100000e+03,   7.09935208e+02],
[ -0.00000000e+00,   7.09935208e+02],
[  1.00400000e+03,   7.07813888e+02],
[  1.00100000e+03,   7.09935208e+02],
[  0.00000000e+00,   7.07813888e+02]])


Finally, FlowDirectorDINF identifies sinks, or local lows.

>>> mg.at_node['flow__sink_flag'].astype(int)
array([1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1])


The flow directors also have the ability to return the flow receiver nodes through a function called direct_flow()

>>> fd = FlowDirectorDINF(mg, 'topographic__elevation')
>>> fd.run_one_step()
array([[ 0, -1],
[ 1, -1],
[ 2, -1],
[ 3, -1],
[ 4, -1],
[ 0, -1],
[ 5,  1],
[ 7, -1],
[ 8, -1],
[ 5, -1],
[ 5, -1],
[11, -1],
[12, -1],
[13, -1],
[14, -1],
[15, -1]])
>>> proportions # doctest: +NORMALIZE_WHITESPACE
array([[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 0.59033447,  0.40966553],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ],
[ 1.        ,  0.        ]])


For each donor node (represented by each row) the proportions should sum to one.

>>> proportions.sum(axis=1)
array([ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,
1.,  1.,  1.])


References

Required Software Citation(s) Specific to this Component

None Listed

Tarboton, D. (1997). A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resources Research 33(2), 309-319. https://dx.doi.org/10.1029/96wr03137

Parameters: grid (ModelGrid) – A grid. surface (field name at node or array of length node, optional) – The surface to direct flow across, default is field at node: topographic__elevation. partition_method (string, optional) – Method for partitioning flow. Options include ‘slope’ (default) and ‘square_root_of_slope’.
__init__(grid, surface='topographic__elevation')[source]
Parameters: grid (ModelGrid) – A grid. surface (field name at node or array of length node, optional) – The surface to direct flow across, default is field at node: topographic__elevation. partition_method (string, optional) – Method for partitioning flow. Options include ‘slope’ (default) and ‘square_root_of_slope’.
direct_flow()[source]

Find flow directions, save to the model grid, and return receivers.

direct_flow() checks for updated boundary conditions, calculates slopes on links, finds basself._surface_valuesel nodes based on the status at node, calculates flow directions, saves results to the grid, and returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_nodes’]

An alternative to direct_flow() is run_one_step() which does the same things but also returns a at-node array of receiver nodes. This array is stored in the grid at: grid[‘node’][‘flow__receiver_nodes’]

run_one_step()[source]

Find flow directions and save to the model grid.

run_one_step() checks for updated boundary conditions, calculates slopes on links, finds basself._surface_valuesel nodes based on the status at node, calculates flow directions, and saves results to the grid.

An alternative to direct_flow() is direct_flow() which does the same things but also returns the receiver nodes not return values.

updated_boundary_conditions()[source]

Method to update FlowDirectorDINF when boundary conditions change.

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