flow_accumulator.py: Component to accumulate flow and calculate drainage area.
Provides the FlowAccumulator component which accumulates flow and calculates drainage area. FlowAccumulator supports multiple methods for calculating flow direction. Optionally a depression finding component can be specified and flow directing, depression finding, and flow routing can all be accomplished together.
FlowAccumulator
(grid, surface='topographic__elevation', flow_director='FlowDirectorSteepest', runoff_rate=None, depression_finder=None, **kwargs)[source]¶Bases: landlab.core.model_component.Component
Component to accumulate flow and calculate drainage area.
This is accomplished by first finding flow directions by a userspecified method and then calculating the drainage area and discharge.
Optionally, spatially variable runoff can be set either by the model grid field ‘water__unit_flux_in’ or the input variable runoff_rate*.
Optionally a depression finding component can be specified and flow directing, depression finding, and flow routing can all be accomplished together.
NOTE: The perimeter nodes 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.
FlowAccumulator stores as ModelGrid fields:
 Node array of drainage areas: ‘drainage_area’
 Node array of discharges: ‘surface_water__discharge’
 Node array containing downstreamtoupstream ordered list of node
 IDs: ‘flow__upstream_node_order’
 Node array of all but the first element of the delta data structure:
 flow__data_structure_delta. The first element is always zero.
The FlowDirector component will add additional ModelGrid fields. DirectToOne methods(Steepest/D4 and D8) and DirectToMany(DINF and MFD) use the same model grid field names. Some of these fields will be different shapes if a DirectToOne or a DirectToMany method is used.
The FlowDirectors store the following 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 for RouteToMany methods and has the shape (nnodes x max number of receivers).
 Node array of flow proportions: ‘flow__receiver_proportions’. This
 array is 2D, for RouteToMany methods and has the shape (nnodes x max number of receivers).
 Node array of links carrying flow: ‘flow__link_to_receiver_node’.
 This array is 2D for RouteToMany methods and has the shape (nnodes x max number of receivers).
 Node array of downhill slopes from each receiver:
 ‘topographic__steepest_slope’ This array is 2D for RouteToMany methods and has the shape (nnodes x max number of receivers).
 Boolean node array of all local lows: ‘flow__sink_flag’
 Link array identifing if flow goes with (1) or against (1) the link direction: ‘flow__link_direction’
The primary method of this class is run_one_step
.
run_one_step takes the optional argument update_flow_director (default is True) that determines if the flow_director is rerun before flow is accumulated.
Parameters: 


Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> mg = RasterModelGrid((3,3))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
The FlowAccumulator component accumulates flow and calculates drainage using all of the different methods for directing flow in Landlab. These include steepest descent (also known as D4 for the case of a raster grid) and D8 (on raster grids only). The method for flow director can be specified as a string (e.g., ‘D8’ or ‘FlowDirectorD8’), as an uninstantiated FlowDirector component or as an instantiated FlowDirector component.
The default method is to use FlowDirectorSteepest.
First let’s look at the three ways to instantiate a FlowAccumulator. The following four methods are all equivalent. First, we can pass the entire name of a flow director as a string to the argument flow_director:
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director='FlowDirectorSteepest'
... )
Second, we can pass just the method name as a string to the argument flow_director:
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director='Steepest'
... )
Third, we can import a FlowDirector component from Landlab and pass it to flow_director:
>>> from landlab.components import FlowDirectorSteepest
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
Finally, we can instantiate a FlowDirector component and pass this instantiated version to flow_director. You might want to do this if you used a FlowDirector in order to set up something before starting a time loop and then want to use the same flow director within the loop.
>>> fd = FlowDirectorSteepest(mg, 'topographic__elevation')
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
Now let’s look at what FlowAccumulator does. Even before we run FlowAccumulator it has the property surface_values that stores the values of the surface over which flow is directed and accumulated.
>>> fa.surface_values
array([ 0., 1., 2., 1., 2., 3., 2., 3., 4.])
Now let’s make a more complicated elevation grid for the next examples.
>>> mg = RasterModelGrid((5, 4))
>>> topographic__elevation = 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("topographic__elevation", topographic__elevation, at="node")
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
>>> fa.run_one_step()
>>> mg.at_node['flow__receiver_node'] # doctest: +NORMALIZE_WHITESPACE
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'] # doctest: +NORMALIZE_WHITESPACE
array([ 0., 1., 5., 0.,
0., 1., 5., 0.,
0., 1., 4., 0.,
0., 1., 2., 0.,
0., 0., 0., 0.])
Now let’s change the cell area (100.) and the runoff rates:
>>> mg = RasterModelGrid((5, 4), xy_spacing=(10., 10))
Put the data back into the new grid.
>>> _ = mg.add_field("topographic__elevation", topographic__elevation, at="node")
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
>>> runoff_rate = np.arange(mg.number_of_nodes, dtype=float)
>>> rnff = mg.add_field("water__unit_flux_in", runoff_rate, at="node", clobber=True)
>>> fa.run_one_step()
>>> mg.at_node['surface_water__discharge'] # doctest: +NORMALIZE_WHITESPACE
array([ 0., 500., 5200., 0.,
0., 500., 5200., 0.,
0., 900., 4600., 0.,
0., 1300., 2700., 0.,
0., 0., 0., 0.])
The flow accumulator will happily work with a negative runoff rate, which could be used to allow, e.g., transmission losses:
>>> runoff_rate.fill(1.)
>>> fa.run_one_step()
>>> mg.at_node['surface_water__discharge']
array([ 0., 100., 500., 0.,
0., 100., 500., 0.,
0., 100., 400., 0.,
0., 100., 200., 0.,
0., 0., 0., 0.])
>>> runoff_rate[:8] = 0.5
>>> fa.run_one_step()
>>> mg.at_node['surface_water__discharge']
array([ 0., 0., 350., 0.,
0., 0., 350., 0.,
0., 100., 400., 0.,
0., 100., 200., 0.,
0., 0., 0., 0.])
The drainage area array is unaffected, as you would expect:
>>> mg.at_node['drainage_area']
array([ 0., 100., 500., 0.,
0., 100., 500., 0.,
0., 100., 400., 0.,
0., 100., 200., 0.,
0., 0., 0., 0.])
The FlowAccumulator component will work for both raster grids and irregular grids. For the example we will use a Hexagonal Model Grid, a special type of Voroni Grid that has regularly spaced hexagonal cells.
>>> from landlab import HexModelGrid
>>> hmg = HexModelGrid((5, 3), xy_of_lower_left=(1., 0.))
>>> _ = hmg.add_field(
... "topographic__elevation",
... hmg.node_x + np.round(hmg.node_y),
... at="node",
... )
>>> fa = FlowAccumulator(
... hmg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
>>> fa.surface_values
array([ 0. , 1. , 2. ,
0.5, 1.5, 2.5, 3.5,
1. , 2. , 3. , 4. , 5. ,
2.5, 3.5, 4.5, 5.5,
3. , 4. , 5. ])
If the FlowDirector you want to use takes keyword arguments and you want to specify it using a string or uninstantiated FlowDirector class, include those keyword arguments when you create FlowAccumulator.
For example, in the case of a raster grid, FlowDirectorMFD can use only orthogonal links, or it can use both orthogonal and diagonal links.
>>> mg = RasterModelGrid((5, 5))
>>> topographic__elevation = mg.node_y+mg.node_x
>>> _ = mg.add_field("topographic__elevation", topographic__elevation, at="node")
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director='MFD',
... diagonals = True
... )
>>> fa.run_one_step()
>>> mg.at_node['flow__receiver_node'] # doctest: +NORMALIZE_WHITESPACE
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],
[ 4, 1, 1, 1, 1, 1, 1, 1],
[ 5, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 5, 1, 1, 1, 0, 1],
[1, 1, 6, 2, 1, 1, 1, 1],
[1, 1, 7, 3, 1, 1, 2, 1],
[ 9, 1, 1, 1, 1, 1, 1, 1],
[10, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 10, 6, 1, 1, 5, 1],
[1, 1, 11, 7, 1, 1, 6, 1],
[1, 1, 12, 8, 1, 1, 7, 1],
[14, 1, 1, 1, 1, 1, 1, 1],
[15, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 15, 11, 1, 1, 10, 1],
[1, 1, 16, 12, 1, 1, 11, 1],
[1, 1, 17, 13, 1, 1, 12, 1],
[19, 1, 1, 1, 1, 1, 1, 1],
[20, 1, 1, 1, 1, 1, 1, 1],
[21, 1, 1, 1, 1, 1, 1, 1],
[22, 1, 1, 1, 1, 1, 1, 1],
[23, 1, 1, 1, 1, 1, 1, 1],
[24, 1, 1, 1, 1, 1, 1, 1]])
>>> mg.at_node['drainage_area'].round(4) # doctest: +NORMALIZE_WHITESPACE
array([ 1.4117, 2.065 , 1.3254, 0.4038, 0. ,
2.065 , 3.4081, 2.5754, 1.3787, 0. ,
1.3254, 2.5754, 2.1716, 1.2929, 0. ,
0.4038, 1.3787, 1.2929, 1. , 0. ,
0. , 0. , 0. , 0. , 0. ])
It may seem odd that there are no round numbers in the drainage area field. This is because flow is directed to all downhill boundary nodes and partitioned based on slope.
To check that flow is conserved, sum along all boundary nodes.
>>> round(sum(mg.at_node['drainage_area'][mg.boundary_nodes]), 4)
9.0
This should be the same as the number of core nodes — as boundary nodes in landlab do not have area.
>>> len(mg.core_nodes)
9
Next, let’s set the dx spacing such that each cell has an area of one.
>>> dx=(2./(3.**0.5))**0.5
>>> hmg = HexModelGrid((5, 3), spacing=dx, xy_of_lower_left=(1.0745, 0.))
>>> _ = hmg.add_field(
... "topographic__elevation",
... hmg.node_x**2 + np.round(hmg.node_y)**2,
... at="node",
... )
>>> fa = FlowAccumulator(
... hmg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
>>> fa.run_one_step()
>>> hmg.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])
>>> np.round(hmg.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 change the cell area (100.) and the runoff rates:
>>> hmg = HexModelGrid((5, 3), spacing=dx * 10.0, xy_of_lower_left=(10.745, 0.))
Put the data back into the new grid.
>>> _ = hmg.add_field(
... "topographic__elevation",
... hmg.node_x ** 2 + np.round(hmg.node_y) ** 2,
... at="node",
... )
>>> fa = FlowAccumulator(
... hmg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
>>> fa.run_one_step()
>>> np.round(hmg.at_node['surface_water__discharge'])
array([ 500., 0., 0.,
200., 500., 200., 0.,
0., 200., 200., 100., 0.,
0., 100., 100., 0.,
0., 0., 0.])
Next, let’s see what happens to a raster grid when there is a depression.
>>> mg = RasterModelGrid((7, 7), xy_spacing=0.5)
>>> z = mg.add_field("topographic__elevation", mg.node_x.copy(), at="node")
>>> z += 0.01 * mg.node_y
>>> mg.at_node['topographic__elevation'].reshape(mg.shape)[2:5, 2:5] *= 0.1
>>> mg.set_closed_boundaries_at_grid_edges(True, True, False, True)
This model grid has a depression in the center.
>>> mg.at_node['topographic__elevation'].reshape(mg.shape)
array([[ 0. , 0.5 , 1. , 1.5 , 2. , 2.5 , 3. ],
[ 0.005 , 0.505 , 1.005 , 1.505 , 2.005 , 2.505 , 3.005 ],
[ 0.01 , 0.51 , 0.101 , 0.151 , 0.201 , 2.51 , 3.01 ],
[ 0.015 , 0.515 , 0.1015, 0.1515, 0.2015, 2.515 , 3.015 ],
[ 0.02 , 0.52 , 0.102 , 0.152 , 0.202 , 2.52 , 3.02 ],
[ 0.025 , 0.525 , 1.025 , 1.525 , 2.025 , 2.525 , 3.025 ],
[ 0.03 , 0.53 , 1.03 , 1.53 , 2.03 , 2.53 , 3.03 ]])
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest
... )
>>> fa.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, 11, 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, 39, 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.5 , 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.5 , 0.25, 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. , 0. ]])
Because of the depression, the flow ‘got stuck’ in the hole in the center of the grid. We can fix this by using a depression finder, such as DepressionFinderAndRouter.
>>> from landlab.components import DepressionFinderAndRouter
We can either run the depression finder separately from the flow accumulator or we can specify the depression finder and router when we instantiate the accumulator and it will run automatically. Similar to specifying the FlowDirector we can provide a depression finder in multiple three ways.
First let’s try running them separately.
>>> df_4 = DepressionFinderAndRouter(mg)
>>> df_4.map_depressions()
>>> mg.at_node['flow__receiver_node'].reshape(mg.shape)
array([[ 0, 1, 2, 3, 4, 5, 6],
[ 7, 7, 16, 17, 18, 11, 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, 39, 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.5 , 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.5 , 0.25, 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. , 0. ]])
Now the flow is routed correctly. The depression finder has properties that including whether there is a lake at the node, which lake is at each node, the outlet node of each lake, and the area of each lake.
>>> df_4.lake_at_node.reshape(mg.shape) # doctest: +NORMALIZE_WHITESPACE
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_4.lake_map.reshape(mg.shape) # doctest: +NORMALIZE_WHITESPACE
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_4.lake_codes # a unique code for each lake present on the grid
array([16])
>>> df_4.lake_outlets # the outlet node of each lake in lake_codes
array([8])
>>> df_4.lake_areas # the area of each lake in lake_codes
array([ 2.25])
Alternatively, we can initialize a flow accumulator with a depression finder specified. Calling run_one_step() will run both the accumulator and the depression finder with one call. For this example, we will pass the class DepressionFinderAndRouter to the parameter depression_finder.
>>> mg = RasterModelGrid((7, 7), xy_spacing=0.5)
>>> z = mg.add_field("topographic__elevation", mg.node_x.copy(), at="node")
>>> z += 0.01 * mg.node_y
>>> mg.at_node['topographic__elevation'].reshape(mg.shape)[2:5, 2:5] *= 0.1
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director='FlowDirectorD8',
... depression_finder=DepressionFinderAndRouter
... )
>>> fa.run_one_step()
This has the same effect of first calling the accumulator and then calling the depression finder.
>>> 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. ]])
The depression finder is stored as part of the flow accumulator, so its properties can be accessed through the depression finder.
>>> fa.depression_finder.lake_at_node.reshape(mg.shape) # doctest: +NORMALIZE_WHITESPACE
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)
>>> fa.depression_finder.lake_map.reshape(mg.shape) # doctest: +NORMALIZE_WHITESPACE
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]])
>>> fa.depression_finder.lake_codes # a unique code for each lake present on the grid
array([16])
>>> fa.depression_finder.lake_outlets # the outlet node of each lake in lake_codes
array([8])
>>> fa.depression_finder.lake_areas # the area of each lake in lake_codes
array([ 2.25])
Finally, note that the DepressionFinderAndRouter takes a keyword argument routing (‘D8’, default; ‘D4’) that sets how connectivity is set between nodes. Similar to our ability to pass keyword arguments to the FlowDirector through FlowAccumulator, we can pass this keyword argument to the DepressionFinderAndRouter component.
>>> fa = FlowAccumulator(
... mg,
... 'topographic__elevation',
... flow_director=FlowDirectorSteepest,
... depression_finder=DepressionFinderAndRouter,
... routing='D4'
... )
FlowAccumulator was designed to work with all types of grids. However,
NetworkModelGrid’s have no cell area. Thus, in order for FlowAccumulator to
this type of grid, an atnode array called cell_area_at_node
must be
present.
>>> from landlab.grid.network import NetworkModelGrid
>>> y_of_node = (0, 1, 2, 2)
>>> x_of_node = (0, 0, 1, 1)
>>> nodes_at_link = ((1, 0), (2, 1), (3, 1))
>>> nmg = NetworkModelGrid((y_of_node, x_of_node), nodes_at_link)
>>> area = nmg.add_ones("cell_area_at_node", at="node")
>>> z = nmg.add_field(
... "topographic__elevation",
... nmg.x_of_node + nmg.y_of_node,
... at="node",
... )
>>> fa = FlowAccumulator(nmg)
>>> fa.run_one_step()
>>> nmg.at_node['flow__receiver_node']
array([0, 0, 2, 1])
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
Braun, J., Willett, S. (2013). A very efficient O(n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution. Geomorphology 180181(C), 170179. https://dx.doi.org/10.1016/j.geomorph.2012.10.008
Initialize the FlowAccumulator component.
Saves the grid, tests grid type, tests imput types and compatability for the flow_director and depression_finder keyword arguments, tests the argument of runoff_rate, and initializes new fields.
__init__
(grid, surface='topographic__elevation', flow_director='FlowDirectorSteepest', runoff_rate=None, depression_finder=None, **kwargs)[source]¶Initialize the FlowAccumulator component.
Saves the grid, tests grid type, tests imput types and compatability for the flow_director and depression_finder keyword arguments, tests the argument of runoff_rate, and initializes new fields.
accumulate_flow
(update_flow_director=True, update_depression_finder=True)[source]¶Function to make FlowAccumulator calculate drainage area and discharge.
Running accumulate_flow() results in the following to occur:
 Flow directions are updated (unless update_flow_director is set as False). This incluldes checking for updated boundary conditions.
 The depression finder, if present is updated (unless update_depression_finder is set as False).
 Intermediate steps that analyse the drainage network topology and create datastructures for efficient drainage area and discharge calculations.
 Calculation of drainage area and discharge.
 Return of drainage area and discharge.
Parameters: 


Returns: 

depression_finder
¶The DepressionFinder used internally.
flow_director
¶The FlowDirector used internally.
headwater_nodes
()[source]¶Return the headwater nodes.
These are nodes that contribute flow and have no upstream nodes.
Examples
>>> from numpy.testing import assert_array_equal
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> mg = RasterModelGrid((5, 5))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
>>> fa = FlowAccumulator(mg, 'topographic__elevation')
>>> fa.run_one_step()
>>> assert_array_equal(fa.headwater_nodes(), np.array([16, 17, 18]))
link_order_upstream
()[source]¶Return the upstream order of active links.
Examples
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> mg = RasterModelGrid((5, 5))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
>>> fa = FlowAccumulator(mg, 'topographic__elevation')
>>> fa.run_one_step()
>>> fa.link_order_upstream()
array([ 5, 14, 23, 6, 15, 24, 7, 16, 25])
This also works for routetomany methods
>>> mg = RasterModelGrid((5, 5))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, True, False)
>>> _ = mg.add_field(
... "topographic__elevation",
... mg.node_x + mg.node_y,
... at="node",
... )
>>> fa = FlowAccumulator(mg,
... 'topographic__elevation',
... flow_director='MFD')
>>> fa.run_one_step()
>>> fa.link_order_upstream()
array([ 5, 14, 10, 6, 11, 7, 23, 19, 15, 20, 16, 28, 24, 29, 25])
node_drainage_area
¶Return the drainage area.
node_order_upstream
¶Return the upstream node order (drainage stack).
node_water_discharge
¶Return the surface water discharge.
run_one_step
()[source]¶Accumulate flow and save to the model grid.
An alternative to run_one_step() is accumulate_flow() which does the same things but also returns the drainage area and discharge. accumulate_flow() additionally provides the ability to turn off updating the flow director or the depression finder.
surface_values
¶Values of the surface over which flow is accumulated.
lossy_flow_accumulator.py: Component to accumulate flow and calc drainage area, while permitting gain or loss of discharge during flow.
DEJH, late 2018
LossyFlowAccumulator
(grid, surface='topographic__elevation', flow_director='FlowDirectorSteepest', runoff_rate=None, depression_finder=None, loss_function=None, **kwargs)[source]¶Bases: landlab.components.flow_accum.flow_accumulator.FlowAccumulator
Component to calculate drainage area and accumulate flow, while permitting dynamic loss or gain of flow downstream.
This component is closely related to the FlowAccumulator, in that this is accomplished by first finding flow directions by a userspecified method and then calculating the drainage area and discharge. However, this component additionally requires the passing of a function that describes how discharge is lost or gained downstream, f(Qw, nodeID, linkID, grid). See the Examples below to see how this works in practice.
Optionally, spatially variable runoff can be set either by the model grid field ‘water__unit_flux_in’ or the input variable runoff_rate*.
Optionally a depression finding component can be specified and flow directing, depression finding, and flow routing can all be accomplished together. Note that the DepressionFinderAndRouter is not particularly intelligent when running on lossy streams, and in particular, it will reroute flow around pits even when they are in fact not filled due to loss.
NOTE: The perimeter nodes 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.
LossyFlowAccumulator stores as ModelGrid fields:
 Node array of drainage areas: ‘drainage_area’
 Node array of discharges: ‘surface_water__discharge’
 Node array of discharge loss in transit (vol/sec). This is the
 total loss across all of the downstream links: ‘surface_water__discharge_loss’
 Node array containing downstreamtoupstream ordered list of node IDs: ‘flow__upstream_node_order’
 Node array of all but the first element of the delta data structure:
 flow__data_structure_delta. The first element is always zero.
The FlowDirector component will add additional ModelGrid fields; see the FlowAccumulator component for full details. These are:
 Node array of receivers (nodes that receive flow), or ITS OWN ID if
 there is no receiver: ‘flow__receiver_node’
 Node array of flow proportions: ‘flow__receiver_proportions’
 Node array of links carrying flow: ‘flow__link_to_receiver_node’
 Node array of downhill slopes from each receiver:
 ‘topographic__steepest_slope’
 Boolean node array of all local lows: ‘flow__sink_flag’
The primary method of this class is run_one_step
.
Examples
These examples pertain only to the LossyFlowAccumulator. See the main FlowAccumulator documentation for more generic and comprehensive examples.
First, a very simple example. Here’s a 50% loss of discharge every time flow moves along a node:
>>> import numpy as np
>>> from landlab import RasterModelGrid, HexModelGrid
>>> from landlab.components import FlowDirectorSteepest
>>> from landlab.components import DepressionFinderAndRouter
>>> mg = RasterModelGrid((3, 5), xy_spacing=(2, 1))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, False, True)
>>> z = mg.add_field("topographic__elevation", mg.node_x + mg.node_y, at="node")
>>> def mylossfunction(qw):
... return 0.5 * qw
>>> fa = LossyFlowAccumulator(mg, 'topographic__elevation',
... flow_director=FlowDirectorSteepest,
... loss_function=mylossfunction)
>>> fa.run_one_step()
>>> mg.at_node['drainage_area'].reshape(mg.shape)
array([[ 0., 0., 0., 0., 0.],
[ 6., 6., 4., 2., 0.],
[ 0., 0., 0., 0., 0.]])
>>> mg.at_node['surface_water__discharge'].reshape(mg.shape)
array([[ 0. , 0. , 0. , 0. , 0. ],
[ 1.75, 3.5 , 3. , 2. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ]])
>>> mg.at_node['surface_water__discharge_loss'].reshape(mg.shape)
array([[ 0. , 0. , 0. , 0. , 0. ],
[ 0. , 1.75, 1.5 , 1. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ]])
Here we use a spatially distributed field to derive loss terms, and also use a filled, nonraster grid.
>>> dx=(2./(3.**0.5))**0.5 # area to be 100.
>>> hmg = HexModelGrid((5, 3), spacing=dx, xy_of_lower_left=(1.0745, 0.))
>>> z = hmg.add_field(
... "topographic__elevation",
... hmg.node_x**2 + np.round(hmg.node_y)**2,
... at="node",
... )
>>> z[9] = 10. # poke a hole
>>> lossy = hmg.add_zeros('node', 'mylossterm', dtype=float)
>>> lossy[14] = 1. # suppress all flow from node 14
Without loss looks like this:
>>> fa = LossyFlowAccumulator(hmg, 'topographic__elevation',
... flow_director=FlowDirectorSteepest,
... depression_finder=DepressionFinderAndRouter)
>>> fa.run_one_step()
>>> hmg.at_node['flow__receiver_node']
array([ 0, 1, 2,
3, 0, 9, 6,
7, 9, 4, 9, 11,
12, 9, 9, 15,
16, 17, 18])
>>> np.round(hmg.at_node['drainage_area'])
array([ 7., 0., 0.,
0., 7., 1., 0.,
0., 1., 6., 1., 0.,
0., 1., 1., 0.,
0., 0., 0.])
>>> np.round(hmg.at_node['surface_water__discharge'])
array([ 7., 0., 0.,
0., 7., 1., 0.,
0., 1., 6., 1., 0.,
0., 1., 1., 0.,
0., 0., 0.])
With loss looks like this:
>>> def mylossfunction2(Qw, nodeID, linkID, grid):
... return (1.  grid.at_node['mylossterm'][nodeID]) * Qw
>>> fa = LossyFlowAccumulator(hmg, 'topographic__elevation',
... flow_director=FlowDirectorSteepest,
... depression_finder=DepressionFinderAndRouter,
... loss_function=mylossfunction2)
>>> fa.run_one_step()
>>> np.round(hmg.at_node['drainage_area'])
array([ 7., 0., 0.,
0., 7., 1., 0.,
0., 1., 6., 1., 0.,
0., 1., 1., 0.,
0., 0., 0.])
>>> np.round(hmg.at_node['surface_water__discharge'])
array([ 6., 0., 0.,
0., 6., 1., 0.,
0., 1., 5., 1., 0.,
0., 1., 1., 0.,
0., 0., 0.])
>>> np.allclose(
... hmg.at_node["surface_water__discharge_loss"],
... lossy * hmg.at_node["surface_water__discharge"],
... )
True
(Loss is only happening from the node, 14, that we set it to happen at.)
Finally, note we can use the linkIDs to create flowlengthdependent effects:
>>> from landlab.components import FlowDirectorMFD
>>> mg = RasterModelGrid((4, 6), xy_spacing=(1, 2))
>>> mg.set_closed_boundaries_at_grid_edges(True, True, False, True)
>>> z = mg.add_field("topographic__elevation", 2.0 * mg.node_x, at="node")
>>> z[9] = 8.
>>> z[16] = 6.5 # force the first node sideways
>>> L = mg.add_zeros('node', 'spatialloss')
>>> mg.at_node['spatialloss'][9] = 1.
>>> mg.at_node['spatialloss'][13] = 1.
>>> def fancyloss(Qw, nodeID, linkID, grid):
... # now a true transmission loss:
... Lt = (1.0  1.0 / grid.length_of_link[linkID] ** 2)
... Lsp = grid.at_node["spatialloss"][nodeID]
... return Qw * (1.0  Lt) * (1.0  Lsp)
>>> fa = LossyFlowAccumulator(
... mg,
... "topographic__elevation",
... flow_director=FlowDirectorMFD,
... loss_function=fancyloss,
... )
>>> fa.run_one_step()
>>> mg.at_node['drainage_area'].reshape(mg.shape)
array([[ 0. , 0. , 0. , 0. , 0. , 0. ],
[ 5.6, 5.6, 3.6, 2. , 2. , 0. ],
[ 10.4, 10.4, 8.4, 6.4, 4. , 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. ]])
>>> mg.at_node['surface_water__discharge'].reshape(mg.shape)
array([[ 0. , 0. , 0. , 0. , 0. , 0. ],
[ 4. , 4. , 2. , 2. , 2. , 0. ],
[ 0. , 8.5, 6.5, 4.5, 2.5, 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. ]])
References
Required Software Citation(s) Specific to this Component
None Listed
Additional References
Braun, J., Willett, S. (2013). A very efficient O(n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution. Geomorphology 180181(C), 170179. https://dx.doi.org/10.1016/j.geomorph.2012.10.008
Initialize the FlowAccumulator component.
Saves the grid, tests grid type, tests imput types and compatability for the flow_director and depression_finder keyword arguments, tests the argument of runoff_rate, and initializes new fields.
Parameters: 


__init__
(grid, surface='topographic__elevation', flow_director='FlowDirectorSteepest', runoff_rate=None, depression_finder=None, loss_function=None, **kwargs)[source]¶Initialize the FlowAccumulator component.
Saves the grid, tests grid type, tests imput types and compatability for the flow_director and depression_finder keyword arguments, tests the argument of runoff_rate, and initializes new fields.
Parameters: 


flow_accum_bw.py: Implementation of the Braun & Willet (2012) stack alorithm.
Implementation of Braun & Willett (2012) algorithm for calculating drainage area and (optionally) water discharge. Assumes each node has only one downstream receiver. If water discharge is calculated, the result assumes steady flow (that is, hydrologic equilibrium).
The main public function is:
a, q, s = flow_accumulation(r)
which takes an array of receivernode IDs, r (the nodes that “receive” the flow from a each node; this array would be returned by the flow_routing component’s calc_flowdirs() method). It returns Numpy arrays with the drainage area (a) and discharge (q) at each node, along with an array (s) that contains the IDs of the nodes in downstreamtoupstream order.
If you simply want the ordered list by itself, use:
s = make_ordered_node_array(r)
Created: GT Nov 2013
find_drainage_area_and_discharge
(s, r, node_cell_area=1.0, runoff=1.0, boundary_nodes=None)[source]¶Calculate the drainage area and water discharge at each node.
Parameters: 


Returns:  drainage area and discharge 
Return type:  tuple of ndarray 
Notes
Examples
>>> import numpy as np
>>> from landlab.components.flow_accum import (
... find_drainage_area_and_discharge)
>>> r = np.array([2, 5, 2, 7, 5, 5, 6, 5, 7, 8])1
>>> s = np.array([4, 1, 0, 2, 5, 6, 3, 8, 7, 9])
>>> a, q = find_drainage_area_and_discharge(s, r)
>>> a
array([ 1., 3., 1., 1., 10., 4., 3., 2., 1., 1.])
>>> q
array([ 1., 3., 1., 1., 10., 4., 3., 2., 1., 1.])
find_drainage_area_and_discharge_lossy
(s, r, link_to_receiver, loss_function, grid, node_cell_area=1.0, runoff=1.0, boundary_nodes=None)[source]¶Calculate the drainage area and water discharge at each node, permitting discharge to fall (or gain) as it moves downstream according to some function. Note that only transmission creates loss, so water sourced locally within a cell is always retained. The loss on each link is recorded in the ‘surface_water__discharge_loss’ link field on the grid; ensure this exists before running the function.
Parameters: 


Returns:  drainage area and discharge 
Return type:  tuple of ndarray 
Notes
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components.flow_accum import (
... find_drainage_area_and_discharge)
>>> r = np.array([2, 5, 2, 7, 5, 5, 6, 5, 7, 8])1
>>> s = np.array([4, 1, 0, 2, 5, 6, 3, 8, 7, 9])
>>> l = np.ones(10, dtype=int) # dummy
>>> nodes_wo_outlet = np.array([0, 1, 2, 3, 5, 6, 7, 8, 9])
>>> def lossfunc(Qw, dummyn, dummyl, dummygrid):
... return 0.5 * Qw
>>> mg = RasterModelGrid((3, 4)) # some grid big enough to make go
>>> _ = mg.add_zeros('node', 'surface_water__discharge_loss', dtype=float)
>>> a, q = find_drainage_area_and_discharge_lossy(s, r, l, lossfunc, mg)
>>> a
array([ 1., 3., 1., 1., 10., 4., 3., 2., 1., 1.])
>>> q
array([ 1. , 2. , 1. , 1. , 3.75, 2. , 2. , 1.5 , 1. , 1. ])
>>> np.allclose(
... mg.at_node['surface_water__discharge_loss'][nodes_wo_outlet],
... 0.5*q[nodes_wo_outlet])
True
>>> np.isclose(mg.at_node['surface_water__discharge_loss'][4], 0.)
True
>>> lossfield = mg.add_ones("loss_field", at="node", dtype=float)
>>> lossfield *= 0.5
>>> def lossfunc2(Qw, nodeID, dummyl, grid):
... return grid.at_node['loss_field'][nodeID] * Qw
>>> a, q = find_drainage_area_and_discharge_lossy(s, r, l, lossfunc2, mg)
>>> a
array([ 1., 3., 1., 1., 10., 4., 3., 2., 1., 1.])
>>> q
array([ 1. , 2. , 1. , 1. , 3.75, 2. , 2. , 1.5 , 1. , 1. ])
>>> np.allclose(
... mg.at_node['surface_water__discharge_loss'][nodes_wo_outlet],
... lossfield[nodes_wo_outlet] * q[nodes_wo_outlet])
True
>>> def lossfunc3(Qw, nodeID, dummyl, dummygrid):
... return Qw  100. # a huge loss
>>> a, q = find_drainage_area_and_discharge_lossy(s, r, l, lossfunc3, mg)
>>> a
array([ 1., 3., 1., 1., 10., 4., 3., 2., 1., 1.])
>>> q
array([ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
flow_accumulation
(receiver_nodes, node_cell_area=1.0, runoff_rate=1.0, boundary_nodes=None)[source]¶Calculate drainage area and (steady) discharge.
Calculates and returns the drainage area and (steady) discharge at each node, along with a downstreamtoupstream ordered list (array) of node IDs.
Examples
>>> import numpy as np
>>> from landlab.components.flow_accum import flow_accumulation
>>> r = np.array([2, 5, 2, 7, 5, 5, 6, 5, 7, 8])1
>>> a, q, s = flow_accumulation(r)
>>> a
array([ 1., 3., 1., 1., 10., 4., 3., 2., 1., 1.])
>>> q
array([ 1., 3., 1., 1., 10., 4., 3., 2., 1., 1.])
>>> s
array([4, 1, 0, 2, 5, 6, 3, 8, 7, 9])
make_ordered_node_array
(receiver_nodes, nd=None, delta=None, D=None)[source]¶Create an array of node IDs that is arranged in order from.
Creates and returns an array of node IDs that is arranged in order from downstream to upstream.
The lack of a leading underscore is meant to signal that this operation could be useful outside of this module!
Examples
>>> import numpy as np
>>> from landlab.components.flow_accum import make_ordered_node_array
>>> r = np.array([2, 5, 2, 7, 5, 5, 6, 5, 7, 8])1
>>> s = make_ordered_node_array(r)
>>> s
array([4, 1, 0, 2, 5, 6, 3, 8, 7, 9])
Short description.
flow_accum_to_n.py: Implementation a routetomultiple drainage stack alorithm.
Algorithm for route to multiple (N) flow accumulation. Inspiration for data structures and attempting O(n) efficiency taken from Braun and Willet(2013).
Algorithm constructs drainage area and (optionally) water discharge. Can handle the case in which each node has more than one downstream receiver.
Computationally, for a grid of the same size this algorithm will take about
 1.5 x (avg number of downstream nodes per cell)
 x (duration of flow_accum_bw for same grid using routetoone method)
So under routetoone direction schemes, using the Braun and Willet method is recommended.
If water discharge is calculated, the result assumes steady flow (that is, hydrologic equilibrium).
The main public function is:
a, q, s = flow_accumulation_to_n(r, p)
which takes the following inputs:
r, an (np, q) array of receivernode IDs, where np is the total number of nodes and q is the maximum number of receivers any node in the grid has. This array would be returned by the flow_routing component.
p, an (np, q) array that identifies the proportion of flow going to each receiver. For each q elements along the np axis, sum(p(i, :)) must equal 1. This array would be returned by the flow_routing component.
It returns Numpy arrays with the drainage area (a) and discharge (q) at each node, along with an array (s) that contains the IDs of the nodes in downstream toupstream order.
If you simply want the ordered list by itself, use:
s = make_ordered_node_array_to_n(r, p, b)
Created: KRB Oct 2016 (modified from flow_accumu_bw)
find_drainage_area_and_discharge_to_n
(s, r, p, node_cell_area=1.0, runoff=1.0, boundary_nodes=None)[source]¶Calculate the drainage area and water discharge at each node.
Parameters: 


Returns:  drainage area and discharge 
Return type:  tuple of ndarray 
Notes
Examples
>>> import numpy as np
>>> from landlab.components.flow_accum.flow_accum_to_n import(
... find_drainage_area_and_discharge_to_n)
>>> r = np.array([[ 1, 2],
... [ 4, 5],
... [ 1, 5],
... [ 6, 2],
... [ 4, 1],
... [ 4, 1],
... [ 5, 7],
... [ 4, 5],
... [ 6, 7],
... [ 7, 8]])
>>> p = np.array([[ 0.6, 0.4 ],
... [ 0.85, 0.15],
... [ 0.65, 0.35],
... [ 0.9, 0.1 ],
... [ 1., 0. ],
... [ 1., 0. ],
... [ 0.75, 0.25],
... [ 0.55, 0.45],
... [ 0.8, 0.2 ],
... [ 0.95, 0.05]])
>>> s = np.array([4, 5, 1, 7, 2, 6, 0, 8, 3, 9])
>>> a, q = find_drainage_area_and_discharge_to_n(s, r, p)
>>> a.round(4)
array([ 1. , 2.575 , 1.5 , 1. , 10. , 5.2465,
2.74 , 2.845 , 1.05 , 1. ])
>>> q.round(4)
array([ 1. , 2.575 , 1.5 , 1. , 10. , 5.2465,
2.74 , 2.845 , 1.05 , 1. ])
find_drainage_area_and_discharge_to_n_lossy
(s, r, link_to_receiver, p, loss_function, grid, node_cell_area=1.0, runoff=1.0, boundary_nodes=None)[source]¶Calculate the drainage area and water discharge at each node, permitting discharge to fall (or gain) as it moves downstream according to some function. Note that only transmission creates loss, so water sourced locally within a cell is always retained. The loss on each link is recorded in the ‘surface_water__discharge_loss’ link field on the grid; ensure this exists before running the function.
Parameters: 


Returns:  drainage area and discharge 
Return type:  tuple of ndarray 
Notes
Examples
>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components.flow_accum.flow_accum_to_n import(
... find_drainage_area_and_discharge_to_n_lossy)
>>> r = np.array([[ 1, 2],
... [ 3, 1],
... [ 3, 1],
... [ 3, 1]])
>>> p = np.array([[ 0.5, 0.5],
... [ 1. , 0. ],
... [ 0.2, 0.8],
... [ 1. , 0. ]])
>>> s = np.array([ 3, 1, 2, 0])
>>> l = np.ones_like(r, dtype=int) # dummy
Make here a grid that contains (too many!) links holding values for loss. We’re only going to use the first 4 links, but illustrates the use of the grid for link input.
>>> mg = RasterModelGrid((3, 3))
>>> _ = mg.add_zeros('node', 'surface_water__discharge_loss', dtype=float)
>>> lossy = mg.add_ones("lossy", at="link", dtype=float)
>>> lossy *= 0.5
>>> def lossfunc(Qw, dummyn, linkID, grid):
... return grid.at_link['lossy'][linkID] * Qw
>>> a, q = find_drainage_area_and_discharge_to_n_lossy(
... s, r, l, p, lossfunc, mg)
>>> a
array([ 1. , 2.7, 1.5, 4. ])
>>> q
array([ 1. , 1.75, 1.25, 2. ])
>>> np.allclose(mg.at_node['surface_water__discharge_loss'][:3], 0.5*q[:3])
True
Note by definition no loss is occuring at the outlet node, as there are no nodes downstream.
Final example of total transmission loss:
>>> def lossfunc(Qw, dummyn, dummyl, dummygrid):
... return Qw  100. # huge loss
>>> a, q = find_drainage_area_and_discharge_to_n_lossy(
... s, r, l, p, lossfunc, mg)
>>> a
array([ 1. , 2.7, 1.5, 4. ])
>>> q
array([ 1., 1., 1., 1.])
flow_accumulation_to_n
(receiver_nodes, receiver_proportions, node_cell_area=1.0, runoff_rate=1.0, boundary_nodes=None)[source]¶Calculate drainage area and (steady) discharge.
Calculates and returns the drainage area and (steady) discharge at each node, along with a downstreamtoupstream ordered list (array) of node IDs.
Examples
>>> import numpy as np
>>> from landlab.components.flow_accum.flow_accum_to_n import(
... flow_accumulation_to_n)
>>> r = np.array([[ 1, 2],
... [ 4, 5],
... [ 1, 5],
... [ 6, 2],
... [ 4, 1],
... [ 4, 1],
... [ 5, 7],
... [ 4, 5],
... [ 6, 7],
... [ 7, 8]])
>>> p = np.array([[ 0.6, 0.4 ],
... [ 0.85, 0.15],
... [ 0.65, 0.35],
... [ 0.9, 0.1 ],
... [ 1., 0. ],
... [ 1., 0. ],
... [ 0.75, 0.25],
... [ 0.55, 0.45],
... [ 0.8, 0.2 ],
... [ 0.95, 0.05]])
>>> a, q, s = flow_accumulation_to_n(r, p)
>>> a.round(4)
array([ 1. , 2.575 , 1.5 , 1. , 10. , 5.2465,
2.74 , 2.845 , 1.05 , 1. ])
>>> q.round(4)
array([ 1. , 2.575 , 1.5 , 1. , 10. , 5.2465,
2.74 , 2.845 , 1.05 , 1. ])
>>> s[0] == 4
True
>>> s[1] == 5
True
>>> s[9] == 9
True
>>> len(set([1, 7])set(s[2:4]))
0
>>> len(set([2, 6])set(s[4:6]))
0
>>> len(set([0, 3, 8])set(s[6:9]))
0
make_ordered_node_array_to_n
(receiver_nodes, receiver_proportion, nd=None, delta=None, D=None)[source]¶Create an array of node IDs.
Creates and returns an array of node IDs that is arranged in order from downstream to upstream.
The lack of a leading underscore is meant to signal that this operation could be useful outside of this module!
Examples
>>> import numpy as np
>>> from landlab.components.flow_accum.flow_accum_to_n import(
... make_ordered_node_array_to_n)
>>> r = np.array([[ 1, 2],
... [ 4, 5],
... [ 1, 5],
... [ 6, 2],
... [ 4, 1],
... [ 4, 1],
... [ 5, 7],
... [ 4, 5],
... [ 6, 7],
... [ 7, 8]])
>>> p = np.array([[ 0.6, 0.4 ],
... [ 0.85, 0.15],
... [ 0.65, 0.35],
... [ 0.9, 0.1 ],
... [ 1., 0. ],
... [ 1., 0. ],
... [ 0.75, 0.25],
... [ 0.55, 0.45],
... [ 0.8, 0.2 ],
... [ 0.95, 0.05]])
>>> s = make_ordered_node_array_to_n(r, p)
>>> s[0] == 4
True
>>> s[1] == 5
True
>>> s[9] == 9
True
>>> len(set([1, 7])set(s[2:4]))
0
>>> len(set([2, 6])set(s[4:6]))
0
>>> len(set([0, 3, 8])set(s[6:9]))
0