Calculate singlepath (steepest direction) flow directions.
Given a ModelGrid, calculates singlepath (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).
FlowRouter
(self, grid, method='D8', runoff_rate=None, **kwds)[source]¶Bases: landlab.components.flow_accum.flow_accumulator.FlowAccumulator
Singlepath (steepest direction) flow routing.
This class implements singlepath (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()
.
Initialize FlowDirector.
Parameters:  grid : ModelGrid
method : {‘D8’, ‘D4’}, optional
runoff_rate : float, optional (m/time)


link_to_flow_receiving_node
¶node_drainage_area
¶node_order_upstream
¶node_receiving_flow
¶node_steepest_slope
¶node_water_discharge
¶route_flow
(**kwds)[source]¶Route surfacewater flow over a landscape.
Routes surfacewater 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:
Returns:  ModelGrid


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 surfacewater flow over a landscape.
Routes surfacewater 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.])
The default behavior of FlowRouter is to use the D8 method. Next we will examine the alternative case of the D4 method that does not consider diagonal links bewtween nodes.
>>> 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, method = 'D4')
>>> fr.run_one_step()
>>> mg.at_node['flow__receiver_node']
array([ 0, 1, 2, 3,
4, 1, 2, 7,
8, 10, 6, 11,
12, 14, 10, 15,
16, 17, 18, 19])
>>> mg.at_node['drainage_area']
array([ 0., 1., 5., 0.,
0., 1., 5., 0.,
0., 1., 4., 0.,
0., 1., 2., 0.,
0., 0., 0., 0.])
The flow router can also work on irregular grids. For the example we will use a Hexagonal Model Grid, a special type of Voroni Grid that has regularly spaced hexagonal cells. We will also set the dx spacing such that each cell has an area of one.
>>> from landlab import HexModelGrid
>>> dx=(2./(3.**0.5))**0.5
>>> mg = HexModelGrid(5,3, dx)
>>> _ = mg.add_field('topographic__elevation', mg.node_x + np.round(mg.node_y), at = 'node')
>>> fr = FlowRouter(mg)
>>> fr.run_one_step()
>>> mg.at_node['flow__receiver_node']
array([ 0, 1, 2,
3, 0, 1, 6,
7, 3, 4, 5, 11,
12, 8, 9, 15,
16, 17, 18])
>>> mg.at_node['drainage_area']
array([ 3., 2., 0.,
2., 3., 2., 0.,
0., 2., 2., 1., 0.,
0., 1., 1., 0.,
0., 0., 0.])
Now let’s return to the first example and 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.])
Find depressions on a topographic surface.
Code author: gtucker, DEJH (Flow routing)
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 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.
Create a DepressionFinderAndRouter.
Constructor assigns a copy of the grid, sets the current time, and calls the initialize method.
Parameters:  grid : RasterModelGrid
routing : ‘D8’ or ‘D4’ (optional)


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])
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


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


Returns:  int

is_valid_outlet
(the_node)[source]¶Check if a node is a valid outlet for the depression.
Parameters:  the_node : int
nodes_this_depression : array_like of int


Returns:  boolean

lake_areas
¶A nlakeslong 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 (nonconsecutive) 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 nlakeslong 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
reroute_flow : bool, optional


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
nodes_this_depression : array_like of int


Returns:  boolean

number_of_lakes
¶Return the number of individual lakes.
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., 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(
... [ 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
method : {‘D8’, ‘D4’}, optional
flow_equation : {‘default’, ‘Manning’, ‘Chezy’}, optional
Chezys_C : float (optional)
Mannings_n : float (optional)


discharges_at_links
¶Return the discharges at links.
Note that if diagonal routing, this will return number_of_d8. Otherwise, it will be number_of_links.