Utilities and Decorators#

Additional utilities are located within the core submodule:

Core Landlab Classes

Module contents#

class StablePriorityQueue[source]

Bases: object

Implements a stable priority queue, that tracks insertion order; i.e., this is used to break ties.

See https://docs.python.org/2/library/heapq.html#priority-queue-implementation-notes & https://www.sciencedirect.com/science/article/pii/S0098300413001337

Examples

>>> q = StablePriorityQueue()
>>> q.add_task('b', priority=2)
>>> q.add_task('a', priority=1)
>>> q.add_task(0, priority=0)
>>> q.add_task('c', priority=2)
>>> q.remove_task(0)
>>> q.pop_task()
'a'
>>> q.peek_at_task()
'b'
>>> np.all(q.tasks_currently_in_queue() == np.array(['b', 'c']))
True
>>> q.pop_task()
'b'
>>> np.all(q.tasks_ever_in_queue() == np.array(['b', 'a', '0', 'c']))
True

If only ints or floats are loaded into the array, the _in_queue methods return arrays with the corresponding data types:

>>> q = StablePriorityQueue()
>>> q.add_task(2, priority=2)
>>> q.add_task(1, priority=1)
>>> np.issubdtype(q.tasks_currently_in_queue().dtype, np.integer)
True
>>> q = StablePriorityQueue()
>>> q.add_task(np.pi)
>>> np.issubdtype(q.tasks_currently_in_queue().dtype, np.floating)
True

Popping from (or peeking at) an empty queue will throw a KeyError:

>>> q = StablePriorityQueue()
>>> try:
...     q.pop_task()
... except KeyError:
...     print('No tasks left')
No tasks left

add_task(task, priority=0)[source]: Add a new task or update the priority of an existing task.

peek_at_task()[source]

Return the lowest priority task without removal.

Raise KeyError if empty.

pop_task()[source]

Remove and return the lowest priority task.

Raise KeyError if empty.

remove_task(task)[source]

Mark an existing task as _REMOVED.

Raise KeyError if not found.

tasks_currently_in_queue()[source]: Return array of nodes currently in the queue.

tasks_ever_in_queue()[source]

Return array of all nodes ever added to this queue object.

Repeats are permitted.

add_halo(data, halo=1, halo_value=None)[source]

Add a halo of no data value to data.

Parameters:

data (array-like) – Array to add the halo to.
halo (int, optional) – The size of the halo.
halo_value (float, optional) – Value to fill the halo with. If not provided, the new data will is not initialized.

Examples

>>> import numpy as np
>>> from landlab.utils import add_halo
>>> data = np.array([[1, 2, 3], [4, 5, 6]])
>>> add_halo(data, halo_value=9)
array([[9, 9, 9, 9, 9], [9, 1, 2, 3, 9], [9, 4, 5, 6, 9], [9, 9, 9, 9, 9]])

calculate_window_statistic(grid, field, func, search_radius, calc_on_closed_nodes=True, **kwargs)[source]

Calculate a statistic using a function within a search window.

Note

This only works on grid nodes (not other grid elements e.g. links) for any ModelGrid type.

This utility outputs an array of length equal to the grid’s number of nodes. Each element of the output array represents the node location in the grid. The value of each element is a function of the nodes within the search window surrounding that node location (see the model grid diagram below).

The grid below contains six columns and five rows with cell spacing set to 10 distance units. This utility iteratively evaluates all nodes in the grid. The diagram shows evaluation of node ID 15 (marked x). If the search radius is set to 20, twice the cell spacing, each node marked with a * is within the search window.

· · * · · ·
· * * * · ·
* * x * * ·
· * * * · ·
· · * · · ·

Increasing the search radius to 25 results in the following search window.

· * * * · ·
* * * * * ·
* * x * * ·
* * * * * ·
· * * * · ·

Decreasing the search radius to 15 results in the following search window.

· · · · · ·
· * * * · ·
· * x * · ·
· * * * · ·
· · · · · ·

The input field can be any field assigned to grid nodes (e.g. “topographic__elevation”) and the input function can be any function that acts on the input field (e.g. “np.min” to find the minimum). The input function may be user defined and may contain any number of inputs, which are input as kwargs.

For example, if the input field is “topographic__elevation” and the input function is np.ptp (peak-to-peak, meaning max minus min value), then the output at node 15 will be the maximum elevation within the search window minus the minimum elevation within the search window (also known as relief). The np.percentile function, however, requires not only the input field, but also an input value to define the “q-th percentile” to be calculated. This second input would be added as a kwarg (e.g. q=90) at the end of the inputs for calculate_window_statistic. Both of these scenarios are shown in the examples below.

Parameters:

grid (ModelGrid) – A Landlab ModelGrid.
field (string) – An existing grid field on which to calculate the statistic of interest. Must exist in grid.
func (function) – The function that calculates the window statistic of field. The first parameter of the function must be the values at nodes within the window, which are used used to calculate the statistic for the node under evaluation. Additional parameters of the function can be passed with kwargs.
search_radius (float) – Radius of window within which the statistic is calculated.
calc_on_closed_nodes (boolean, optional) – Toggle calculation over all nodes including closed nodes (True) or all nodes except closed nodes (False).
kwargs (optional) – Keyword arguments passed to func that are additional to the array of node values within the search window.

Returns:

output – Output array containing the calculated values of the statistic. Same length as input field.

Return type:

ndarray

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.utils import window_statistic

>>> grid = RasterModelGrid((5, 6), xy_spacing=10.0)
>>> grid.set_closed_boundaries_at_grid_edges(False, True, False, True)
>>> z = grid.add_zeros("topographic__elevation", at="node")
>>> z += np.arange(len(z))

Calculate relief using np.ptp function.

>>> relief = calculate_window_statistic(
...     grid, "topographic__elevation", np.ptp, search_radius=15
... )
>>> grid.at_node["topographic__elevation"]
array([ 0.,   1.,   2.,   3.,   4.,   5.,
        6.,   7.,   8.,   9.,  10.,  11.,
       12.,  13.,  14.,  15.,  16.,  17.,
       18.,  19.,  20.,  21.,  22.,  23.,
       24.,  25.,  26.,  27.,  28.,  29.])
>>> relief
array([ 7.,   8.,   8.,   8.,   8.,   7.,
       13.,  14.,  14.,  14.,  14.,  13.,
       13.,  14.,  14.,  14.,  14.,  13.,
       13.,  14.,  14.,  14.,  14.,  13.,
        7.,   8.,   8.,   8.,   8.,   7.])

Calculate relief using np.ptp function excluding closed nodes.

>>> relief = calculate_window_statistic(
...     grid,
...     "topographic__elevation",
...     np.ptp,search_radius=15,
...     calc_on_closed_nodes=False,
... )
>>> grid.at_node["topographic__elevation"]
array([ 0.,   1.,   2.,   3.,   4.,   5.,
        6.,   7.,   8.,   9.,  10.,  11.,
       12.,  13.,  14.,  15.,  16.,  17.,
       18.,  19.,  20.,  21.,  22.,  23.,
       24.,  25.,  26.,  27.,  28.,  29.])
>>> relief
array([ nan,  nan,  nan,  nan,  nan,  nan,
         7.,   8.,   8.,   8.,   8.,   7.,
        13.,  14.,  14.,  14.,  14.,  13.,
         7.,   8.,   8.,   8.,   8.,   7.,
        nan,  nan,  nan,  nan,  nan,  nan])

Calculate 90th percentile using np.percentile function and kwargs.

>>> perc_90 = calculate_window_statistic(
...     grid,
...     "topographic__elevation",
...     np.percentile,search_radius=15,
...     calc_on_closed_nodes=False,
...     q=90
... )
>>> grid.at_node["topographic__elevation"]
array([ 0.,   1.,   2.,   3.,   4.,   5.,
        6.,   7.,   8.,   9.,  10.,  11.,
       12.,  13.,  14.,  15.,  16.,  17.,
       18.,  19.,  20.,  21.,  22.,  23.,
       24.,  25.,  26.,  27.,  28.,  29.])
>>> perc_90
array([ nan,  nan,  nan,  nan,  nan,  nan,
       12.7, 13.5, 14.5, 15.5, 16.5, 16.7,
       18.5, 19.2, 20.2, 21.2, 22.2, 22.5,
       18.7, 19.5, 20.5, 21.5, 22.5, 22.7,
        nan,  nan,  nan,  nan,  nan,  nan])

Calculate relief above 90th percentile elevation using a user-defined function and kwargs.

>>> def max_minus_percentile(elev, q):
...     output = np.max(elev) - np.percentile(elev, q)
...     return output
>>> rel_above_90th_perc = calculate_window_statistic(
...     grid,
...     "topographic__elevation",
...     max_minus_percentile,
...     search_radius=15,
...     calc_on_closed_nodes=False,
...     q=90,
... )
>>> grid.at_node["topographic__elevation"]
array([ 0.,   1.,   2.,   3.,   4.,   5.,
        6.,   7.,   8.,   9.,  10.,  11.,
       12.,  13.,  14.,  15.,  16.,  17.,
       18.,  19.,  20.,  21.,  22.,  23.,
       24.,  25.,  26.,  27.,  28.,  29.])
>>> rel_above_90th_perc
array([ nan,  nan,  nan,  nan,  nan,  nan,
        0.3,  0.5,  0.5,  0.5,  0.5,  0.3,
        0.5,  0.8,  0.8,  0.8,  0.8,  0.5,
        0.3,  0.5,  0.5,  0.5,  0.5,  0.3,
        nan,  nan,  nan,  nan,  nan,  nan])

convert_arc_flow_directions_to_landlab_node_ids(grid, flow_dir_arc)[source]

Convert Arc flow_directions to RasterModelGrid node ids.

This function receives flow directions (D8) from ESRI ArcGIS and converts them to Landlab’s RasterModelGrid node id. ESRI ArcGIS D8 flow directions are either of the eight valid output directions relating to the eight adjacent cells into which flow could travel. The valid output directions are powers of 2 starting from 2^0 (1) in the Eastern neighbor going clockwise to 2^7 (128) at Northeastern neighbor. For more information refer ‘https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/ how-flow-direction-works.htm’

Parameters:

grid (RasterModelGrid) – A grid.
flow_dir_arc (ndarray of int, shape (n_nodes, )) – flow directions derived from ESRII ArcGIS.

Returns:

receiver_nodes – downstream node at each node. Note that this array gives the receiver nodes only for the core nodes. For non-core nodes, a zero is used.

Return type:

ndarray of int, shape (n_nodes, )

count_repeated_values(values)[source]

Count how many times in an array values repeat and where they appear.

Return a list of length n that gives the values and indices of repeated values. The first element of the list will be the values and indices of all values that appear once or the first time repeated values appear. The next element, values that repeat twice or more, and so on. Thus, the length of the returned list will be the maximum number that any value is repeated in x.

Parameters:: values (array_like) – Input array to count repeated values.
Returns:: List of tuples of (repeated_values, indices).
Return type:: list of tuple

Examples

For an array that contains no repeated values, this function just returns a copy of x, and the indices to each element.

>>> import numpy as np
>>> from landlab.utils.count_repeats import count_repeated_values
>>> counts = count_repeated_values(np.array([20, 30, 40], dtype=int))
>>> len(counts)
1
>>> counts[0]
(array([20, 30, 40]), array([0, 1, 2]))

If x contains a repeated value, the first element contains all unique values along with their indices. For repeated values, return indices to their first occurance. The second element contains values and indices to values occuring two or more times.

>>> counts = count_repeated_values(np.array([20, 30, 40, 30, 30],
...     dtype=int))
>>> len(counts)
3
>>> counts[0]
(array([20, 30, 40]), array([0, 1, 2]))
>>> counts[1]
(array([30]), array([3]))
>>> counts[2]
(array([30]), array([4]))

The input array remains unchanged.

>>> x = np.array([20, 30, 30, 40], dtype=int)
>>> counts = count_repeated_values(x)
>>> x
array([20, 30, 30, 40])

find_unique_upstream_hsd_ids_and_fractions(hsd_upstr)[source]

Finds unique entries in hsd_upstr.values()

This function operates on hsd_upstr.values(), that are lists of hsd_ids. Two new Python dictionaries, ‘unique_ids’ and ‘fractions’ are created.

unique_ids.keys() = hsd_upstr.keys() unique_ids.values()[i] = list of unique entries in hsd_upstr.values()[i]

fractions.keys() = hsd_upstr.keys() fractions.values()[i] = (number of entries of each unique_id.values()[i]/ length of hsd_upstr.values()[i]) for each unique_id.values()[i] in the same order.

Note that ‘hsd_upstr’ is the output of track_source(). You can use an alternative input. In that case, please refer to the documentation of track_source() or refer source_tracking_algorithm_user_manual for more information.

Parameters:: hsd_upstr (dictionary) – ‘hsd_upstr’ maps each MD grid node to corresponding contributing upstream HSD ids.
Returns:: (unique_ids, fractions) – Tuple of data. ‘unique_ids’ maps each MD node with all upstream HSD ids without repitition. ‘fractions’ maps each MD node with the fractions of contributions of the corresponding upstream HSD ids in the same order as uniques_ids[node_id].
Return type:: (dictionary, dictionary)

get_core_node_at_node(grid)[source]

Get node ids as numbered by core nodes.

Get the core node ID for each node of a grid. If a node is not a core node, then use -1.

Parameters:: grid (ModelGrid) – A ModelGrid.
Returns:: Ids of each of the grid’s core nodes.
Return type:: ndarray of int

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.utils import get_core_node_at_node

>>> grid = RasterModelGrid((4, 5))
>>> get_core_node_at_node(grid)
array([-1, -1, -1, -1, -1,
       -1,  0,  1,  2, -1,
       -1,  3,  4,  5, -1,
       -1, -1, -1, -1, -1])

>>> grid.status_at_node[13] = grid.BC_NODE_IS_FIXED_VALUE
>>> grid.status_at_node[2] = grid.BC_NODE_IS_CLOSED
>>> get_core_node_at_node(grid)
array([-1, -1, -1, -1, -1,
       -1,  0,  1,  2, -1,
       -1,  3,  4, -1, -1,
       -1, -1, -1, -1, -1])

get_core_node_matrix(grid, value_at_node, coef_at_link=None)[source]

A matrix for core nodes and a right-hand side vector.

Construct and return a matrix for the core nodes, plus a right-hand side vector containing values based on the input array value_at_node. Optionally, coef_at_link provides coefficients for each link (default is 1.0).

Parameters:

grid (RasterModelGrid, HexModelGrid) – A landlab grid.
value_at_node (ndarray) – Values defined at nodes used to construct the right-hand side vector.
coef_at_link (ndarray, optional) – Coefficents at links used to construct the matrix. If not provided, use 1.0.

Examples

>>> from landlab import RasterModelGrid
>>> from landlab.utils import get_core_node_matrix

>>> grid = RasterModelGrid((4, 5))
>>> grid.status_at_node[13] = grid.BC_NODE_IS_FIXED_VALUE
>>> grid.status_at_node[2] = grid.BC_NODE_IS_CLOSED

>>> vals = np.arange(grid.number_of_nodes, dtype=np.double)  # made-up state variable array

>>> mat, rhs = get_core_node_matrix(grid, vals)
>>> mat.toarray()
array([[-4.,  1.,  0.,  1.,  0.],
       [ 1., -3.,  1.,  0.,  1.],
       [ 0.,  1., -4.,  0.,  0.],
       [ 1.,  0.,  0., -4.,  1.],
       [ 0.,  1.,  0.,  1., -4.]])
>>> rhs
array([[ -6.],
       [  0.],
       [-25.],
       [-26.],
       [-30.]])

>>> coefs = np.arange(grid.number_of_links, dtype=np.double)  # coefficient array
>>> mat, rhs = get_core_node_matrix(grid, vals, coef_at_link=coefs)
>>> mat.toarray()
array([[-38.,  10.,   0.,  14.,   0.],
       [ 10., -36.,  11.,   0.,  15.],
       [  0.,  11., -46.,   0.,   0.],
       [ 14.,   0.,   0., -74.,  19.],
       [  0.,  15.,   0.,  19., -78.]])
>>> rhs
array([[ -6.],
       [  0.],
       [-25.],
       [-26.],
       [-30.]])

get_watershed_mask(grid, outlet_id)[source]

Get the watershed of an outlet returned as a boolean array.

Parameters:

grid (RasterModelGrid) – A landlab RasterModelGrid.
outlet_id (integer) – The id of the outlet node.

Returns:

watershed_mask – True elements of this array correspond to nodes with flow that is received by the outlet. The length of the array is equal to the grid number of nodes.

Return type:

boolean ndarray

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> from landlab.utils import get_watershed_mask

>>> rmg = RasterModelGrid((7, 7))
>>> z = np.array([
...     -9999., -9999., -9999., -9999., -9999., -9999., -9999.,
...     -9999.,    26.,     0.,    30.,    32.,    34., -9999.,
...     -9999.,    28.,     1.,    25.,    28.,    32., -9999.,
...     -9999.,    30.,     3.,     3.,    11.,    34., -9999.,
...     -9999.,    32.,    11.,    25.,    18.,    38., -9999.,
...     -9999.,    34.,    32.,    34.,    36.,    40., -9999.,
...     -9999., -9999., -9999., -9999., -9999., -9999., -9999.])
>>> rmg.at_node['topographic__elevation'] = z

Only the bottom boundary is set to open. >>> rmg.set_closed_boundaries_at_grid_edges(True, True, True, False) >>> rmg.set_fixed_value_boundaries_at_grid_edges(False, False, False, True)

Route flow. >>> fr = FlowAccumulator(rmg, flow_director=’D8’) >>> fr.run_one_step()

>>> get_watershed_mask(rmg, 2).reshape(rmg.shape)
array([[False, False,  True, False, False, False, False],
       [False, False,  True, False, False, False, False],
       [False,  True,  True,  True, True,  True,  False],
       [False,  True,  True,  True,  True,  True, False],
       [False,  True,  True,  True,  True,  True, False],
       [False,  True,  True,  True,  True,  True, False],
       [False, False, False, False, False, False, False]], dtype=bool)

get_watershed_masks(grid)[source]

Assign the watershed outlet id to all nodes in the grid.

Parameters:: grid (RasterModelGrid) – A landlab RasterModelGrid.
Returns:: watershed_masks – The length of the array is equal to the grid number of nodes. Values of this array are the watershed ids. The value of a watershed id is the node id of the watershed outlet.
Return type:: integer ndarray

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> from landlab.utils import get_watershed_masks

Create a grid with a node spacing of 200 meter. >>> rmg = RasterModelGrid((7, 7), xy_spacing=200) >>> z = np.array([ … -9999., -9999., -9999., -9999., -9999., -9999., -9999., … -9999., 26., 0., 26., 30., 34., -9999., … -9999., 28., 1., 28., 5., 32., -9999., … -9999., 30., 3., 30., 10., 34., -9999., … -9999., 32., 11., 32., 15., 38., -9999., … -9999., 34., 32., 34., 36., 40., -9999., … -9999., -9999., -9999., -9999., -9999., -9999., -9999.]) >>> rmg.at_node[‘topographic__elevation’] = z >>> rmg.set_closed_boundaries_at_grid_edges(True, True, True, False)

Route flow.

>>> fr = FlowAccumulator(rmg, flow_director='D8')
>>> fr.run_one_step()

Assign mask.

>>> mask = get_watershed_masks(rmg)
>>> mask.reshape(rmg.shape)
array([[ 0,  1,  2,  3,  4,  5,  6],
       [ 7,  1,  2,  3,  4,  5, 13],
       [14,  2,  2,  2, 18, 18, 20],
       [21,  2,  2,  2, 18, 18, 27],
       [28,  2,  2,  2, 18, 18, 34],
       [35,  2,  2,  2, 18, 18, 41],
       [42, 43, 44, 45, 46, 47, 48]])

get_watershed_masks_with_area_threshold(grid, critical_area)[source]

Get masks of all of the watersheds with a minimum drainage area size.

Parameters:

grid (RasterModelGrid) – A landlab RasterModelGrid.
critical_area (integer or float) – The minimum drainage area of the watersheds to identify.

Returns:

watershed_masks – The length of the array is equal to the grid number of nodes. Values of this array are the watershed ids. The value of a watershed id is the node id of the watershed outlet. Nodes with a value of -1 have only downstream nodes with drainage areas below critical_area.

Return type:

integer ndarray

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> from landlab.utils import get_watershed_masks_with_area_threshold

Create a grid with a node spacing of 200 meter.

>>> rmg = RasterModelGrid((7, 7), xy_spacing=200)
>>> z = np.array([
...     -9999., -9999., -9999., -9999., -9999., -9999., -9999.,
...     -9999.,    26.,     0.,    26.,    30.,    34., -9999.,
...     -9999.,    28.,     1.,    28.,     5.,    32., -9999.,
...     -9999.,    30.,     3.,    30.,    10.,    34., -9999.,
...     -9999.,    32.,    11.,    32.,    15.,    38., -9999.,
...     -9999.,    34.,    32.,    34.,    36.,    40., -9999.,
...     -9999., -9999., -9999., -9999., -9999., -9999., -9999.])
>>> rmg.at_node['topographic__elevation'] = z
>>> rmg.set_closed_boundaries_at_grid_edges(True, True, True, False)

Route flow.

>>> fr = FlowAccumulator(rmg, flow_director='D8')
>>> fr.run_one_step()

Get the masks of watersheds greater than or equal to 80,000 square-meters.

>>> critical_area = 80000
>>> mask = get_watershed_masks_with_area_threshold(rmg, critical_area)

Verify that all mask null nodes have a drainage area below critical area.

>>> null_nodes = np.where(mask == -1)[0]
>>> A = rmg.at_node['drainage_area'][null_nodes]
>>> below_critical_area_nodes = A < critical_area
>>> np.all(below_critical_area_nodes)
True

get_watershed_nodes(grid, outlet_id)[source]

Get the watershed of an outlet returned as a list of nodes.

Parameters:

grid (RasterModelGrid) – A landlab RasterModelGrid.
outlet_id (integer) – The id of the outlet node.

Returns:

watershed_nodes – The ids of the nodes that flow to the node with the id, outlet_id.

Return type:

integer ndarray

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> from landlab.utils import get_watershed_nodes

>>> rmg = RasterModelGrid((7, 7))
>>> z = np.array([
...     -9999., -9999., -9999., -9999., -9999., -9999., -9999.,
...     -9999.,    26.,     0.,    30.,    32.,    34., -9999.,
...     -9999.,    28.,     1.,    25.,    28.,    32., -9999.,
...     -9999.,    30.,     3.,     3.,    11.,    34., -9999.,
...     -9999.,    32.,    11.,    25.,    18.,    38., -9999.,
...     -9999.,    34.,    32.,    34.,    36.,    40., -9999.,
...     -9999., -9999., -9999., -9999., -9999., -9999., -9999.])

>>> rmg.at_node['topographic__elevation'] = z
>>> rmg.set_watershed_boundary_condition_outlet_id(2, z,
...                                                nodata_value=-9999.)

Route flow. >>> fr = FlowAccumulator(rmg, flow_director=’D8’) >>> fr.run_one_step()

Get the nodes of two watersheds. >>> mainstem_watershed_nodes = get_watershed_nodes(rmg, 2) >>> tributary_watershed_nodes = get_watershed_nodes(rmg, 24)

Given the watershed boundary conditions, the number of mainstem watershed nodes should be equal to the number of core nodes plus the outlet node. >>> len(mainstem_watershed_nodes) == rmg.number_of_core_nodes + 1 True

get_watershed_outlet(grid, source_node_id)[source]

Get the downstream-most node (the outlet) of the source node.

Parameters:

grid (RasterModelGrid) – A landlab RasterModelGrid.
source_node_id (integer) – The id of the node in which to identify its outlet.

Returns:

outlet_node – The id of the node that is the downstream-most node (the outlet) of the source node.

Return type:

integer

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> from landlab.components import FlowAccumulator
>>> from landlab.utils import get_watershed_outlet

>>> rmg = RasterModelGrid((7, 7))
>>> z = np.array([
...     -9999., -9999., -9999., -9999., -9999., -9999., -9999.,
...     -9999.,    26.,     0.,    30.,    32.,    34., -9999.,
...     -9999.,    28.,     1.,    25.,    28.,    32., -9999.,
...     -9999.,    30.,     3.,     3.,    11.,    34., -9999.,
...     -9999.,    32.,    11.,    25.,    18.,    38., -9999.,
...     -9999.,    34.,    32.,    34.,    36.,    40., -9999.,
...     -9999., -9999., -9999., -9999., -9999., -9999., -9999.])

>>> rmg.at_node['topographic__elevation'] = z
>>> imposed_outlet = 2
>>> rmg.set_watershed_boundary_condition_outlet_id(imposed_outlet, z,
...                                                nodata_value=-9999.)

Route flow. >>> fr = FlowAccumulator(rmg, flow_director=’D8’) >>> fr.run_one_step()

Get the grid watershed outlet. >>> determined_outlet = get_watershed_outlet(rmg, 40) >>> determined_outlet == imposed_outlet True

return_array_at_link(grid, value)[source]

Function to return an array stored at node or of shape (n_nodes,).

This function exists to take advantange of the use_field_name_array_or_value decorator which permits providing the surface as a field name or array.

Parameters:

grid (ModelGrid) –
value (field name, ndarray of shape (n_nodes, ), or single value.) –

Returns:

array

Return type:

ndarray of shape (n_nodes, )

return_array_at_node(grid, value)[source]

Function to return an array stored at node or of shape (n_nodes,).

This function exists to take advantange of the use_field_name_array_or_value decorator which permits providing the surface as a field name or array.

Parameters:

grid (ModelGrid) –
value (field name, ndarray of shape (n_nodes, ), or single value.) –

Returns:

array

Return type:

ndarray of shape (n_nodes, )

track_source(grid, hsd_ids, flow_directions=None)[source]

Track all contributing upstream core nodes for each core node.

This algorithm traverses the grid based on information of flow directions at nodes and at every node identifies all the nodes upstream of a given node. The algorithm creates a dictionary with an entry for each node; a node’s entry in the dictionary will contain a list with the node_ids of all upstream nodes. Thus this method permits identification of the source area contributing to each and every node in the model grid. This function is different from a standard flow accumulation routine in that it not only calculates the amount of flow at each node, but records the IDs of all upstream nodes. However, similar to a standard flow accumulation routine, it produces an at_node array of the amount of flow passing through the node. It also differs from a standard flow accumulation routing in that it permits the mapping of flow inputs from a coarser grid to to a finer model grid.

In its present implementation, the algorithm has not been optimized for efficient time use. Its methods are brute force and it should be expected to be time intensive. It is not recommended to be run frequently in a modeling exercise. Due to its intensive nature, this algorithm may fail with large watersheds (a present, the development team has not derived a maximum stable watershed size).

This function was initially developed to find contributing area of a 30 m grid (MD), where the quantitative data that we were interested in was available in significantly coarser resolution (called Hydrologic Source Domain (HSD)). Therefore, we started working with re-sampled HSD, that is at the same resolution as MD, and represents exactly the same landscape. Alternatively, one can use the node ids of MD (grid.nodes.flatten()) as input for hsd_ids.

For more information, refer Ref 1.

Parameters:

grid (RasterModelGrid) – A grid.
hsd_ids (ndarray of int, shape (n_nodes, )) – array that maps the nodes of the grid to, possibly coarser, Hydrologic Source Domain (HSD) grid ids.
flow_directions (ndarray of int, shape (n_nodes, ), optional.) – downstream node at each node. Alternatively, this data can be provided as a nodal field ‘flow__receiver_node’ on the grid.

Returns:

(hsd_upstr, flow_accum) – ‘hsd_upstr’ maps each grid node to corresponding contributing upstream hsd_ids. hsd_upstr.keys() will return node_ids of the grid. hsd_upstr.values() will return lists of all upstream contributing hsd_ids, including repitions of hsd_ids, at corresponding node_ids. ‘flow_accum’ is an array of the number of upstream contributing nodes at each node.

Return type:

(dictionary, ndarray of shape (n_nodes))