landlab

ChiFinder: Perform chi-index analysis for gridded topography

class ChiFinder(grid, reference_concavity=0.5, min_drainage_area=1000000.0, reference_area=1.0, use_true_dx=False, **kwds)[source]

Bases: landlab.core.model_component.Component

This component calculates chi indices, sensu Perron & Royden, 2013, for a Landlab landscape.

Construction:

ChiFinder(grid, reference_concavity=0.5, min_drainage_area=1.e6,
          reference_area=1., use_true_dx=False)
Parameters:

grid : RasterModelGrid

A landlab RasterModelGrid.

reference_concavity : float

The reference concavity to use in the calculation.

min_drainage_area : float (m**2)

The drainage area down to which to calculate chi.

reference_area : float or None (m**2)

If None, will default to the mean core cell area on the grid. Else, provide a value to use. Essentially becomes a prefactor on the value of chi.

use_true_dx : bool (default False)

If True, integration to give chi is performed using each value of node spacing along the channel (which can lead to a quantization effect, and is not preferred by Taylor & Royden). If False, the mean value of node spacing along the all channels is assumed everywhere.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid, CLOSED_BOUNDARY
>>> from landlab.components import FlowRouter, FastscapeEroder
>>> mg = RasterModelGrid((3, 4), 1.)
>>> for nodes in (mg.nodes_at_right_edge, mg.nodes_at_bottom_edge,
...               mg.nodes_at_top_edge):
...     mg.status_at_node[nodes] = CLOSED_BOUNDARY
>>> _ = mg.add_field('node', 'topographic__elevation', mg.node_x)
>>> fr = FlowRouter(mg)
>>> cf = ChiFinder(mg, min_drainage_area=1., reference_concavity=1.)
>>> fr.run_one_step()
>>> cf.calculate_chi()
>>> mg.at_node['channel__chi_index'].reshape(mg.shape)[1, :]
array([ 0.5,  1. ,  2. ,  0. ])
>>> mg2 = RasterModelGrid((5, 5), 100.)
>>> for nodes in (mg2.nodes_at_right_edge, mg2.nodes_at_bottom_edge,
...               mg2.nodes_at_top_edge):
...     mg2.status_at_node[nodes] = CLOSED_BOUNDARY
>>> _ = mg2.add_zeros('node', 'topographic__elevation')
>>> mg2.at_node['topographic__elevation'][mg2.core_nodes] = mg2.node_x[
...     mg2.core_nodes]/1000.
>>> np.random.seed(0)
>>> mg2.at_node['topographic__elevation'][
...     mg2.core_nodes] += np.random.rand(mg2.number_of_core_nodes)
>>> fr2 = FlowRouter(mg2)
>>> sp2 = FastscapeEroder(mg2, K_sp=0.01)
>>> cf2 = ChiFinder(mg2, min_drainage_area=0., reference_concavity=0.5)
>>> for i in range(10):
...     mg2.at_node['topographic__elevation'][mg2.core_nodes] += 10.
...     fr2.run_one_step()
...     sp2.run_one_step(1000.)
>>> fr2.run_one_step()
>>> cf2.calculate_chi()
>>> mg2.at_node['channel__chi_index'].reshape(
...     mg2.shape)  
array([[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.77219416,  1.54438833,  2.63643578,  2.61419437,  0.        ],
       [ 1.09204746,  2.18409492,  1.52214691,  2.61419437,  0.        ],
       [ 0.44582651,  0.89165302,  1.66384718,  2.75589464,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ]])
>>> cf2.calculate_chi(min_drainage_area=20000., use_true_dx=True,
...                   reference_area=mg2.at_node['drainage_area'].max())
>>> cf2.chi_indices.reshape(mg2.shape)  
array([[   0. ,   0.        ,   0.        ,   0. ,   0. ],
       [   0. , 173.20508076,   0.        ,   0. ,   0. ],
       [   0. ,   0.        , 270.71067812,   0. ,   0. ],
       [   0. , 100.        , 236.60254038,   0. ,   0. ],
       [   0. ,   0.        ,   0.        ,   0. ,   0. ]])
>>> cf2.hillslope_mask.reshape(mg2.shape)
array([[ True,  True,  True,  True,  True],
       [False, False,  True,  True,  True],
       [ True,  True, False,  True,  True],
       [False, False, False,  True,  True],
       [ True,  True,  True,  True,  True]], dtype=bool)
best_fit_chi_elevation_gradient_and_intercept(ch_nodes=None)[source]

Returns least squares best fit for a straight line through a chi plot.

Parameters:

ch_nodes : array of ints or None

Nodes at which to consider chi and elevation values. If None, will use all nodes in grid with area greater than the component min_drainage_area.

Returns:

coeffs : array(gradient, intercept)

A len-2 array containing the m then z0, where z = z0 + m * chi.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid, CLOSED_BOUNDARY
>>> from landlab.components import FlowRouter
>>> mg = RasterModelGrid((3, 4), 1.)
>>> for nodes in (mg.nodes_at_right_edge, mg.nodes_at_bottom_edge,
...               mg.nodes_at_top_edge):
...     mg.status_at_node[nodes] = CLOSED_BOUNDARY
>>> z = mg.add_field('node', 'topographic__elevation',
...                  mg.node_x.copy())
>>> z[4:8] = np.array([0.5, 1., 2., 0.])
>>> fr = FlowRouter(mg)
>>> cf = ChiFinder(mg, min_drainage_area=1., reference_concavity=1.)
>>> fr.run_one_step()
>>> cf.calculate_chi()
>>> mg.at_node['channel__chi_index'].reshape(mg.shape)[1, :]
array([ 0.5,  1. ,  2. ,  0. ])
>>> coeffs = cf.best_fit_chi_elevation_gradient_and_intercept()
>>> np.allclose(np.array([1., 0.]), coeffs)
True
calculate_chi(**kwds)[source]

This is the main method. Call it to calculate local chi indices at all points with drainage areas greater than min_drainage_area.

This “run” method can optionally take the same parameter set as provided at instantiation. If they are provided, they will override the existing values from instantiation.

Chi of any node without a defined value is reported as 0. These nodes are also identified in the mask retrieved with hillslope_mask().

chi_indices

Return the array of channel steepness indices.

Nodes not in the channel receive zeros.

create_chi_plot(channel_heads=None, label_axes=True, symbol='kx', plot_line=False, line_symbol='r-')[source]

Plots a “chi plot” (chi vs elevation for points in channel network).

If channel_heads is provided, only the channel nodes downstream of the provided points (and with area > min_drainage_area) will be plotted.

Parameters:

channel_heads : int, list or array of ints, or None

Node IDs of channel heads to from which plot downstream.

label_axes : bool

If True, labels the axes as “Chi” and “Elevation (m)”.

symbol : str

A matplotlib-style string for the style to use for the points.

plot_line : bool

If True, will plot a linear best fit line through the data cloud.

line_symbol : str

A matplotlib-style string for the style to use for the line, if plot_line.

hillslope_mask

Return a boolean array, False where steepness indices exist.

integrate_chi_avg_dx(valid_upstr_order, chi_integrand, chi_array, mean_dx)[source]

Calculates chi at each channel node by summing chi_integrand.

This method assumes a uniform, mean spacing between nodes. Method is deliberately split out for potential cythonization at a later stage.

Parameters:

valid_upstr_order : array of ints

nodes in the channel network in upstream order.

chi_integrand : array of floats

The value (A0/A)**concavity, in upstream order.

chi_array : array of floats

Array in which to store chi.

mean_dx : float

The mean node spacing in the network.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid, CLOSED_BOUNDARY
>>> from landlab.components import FlowRouter
>>> mg = RasterModelGrid((5, 4), 1.)
>>> for nodes in (mg.nodes_at_right_edge, mg.nodes_at_bottom_edge,
...               mg.nodes_at_top_edge):
...     mg.status_at_node[nodes] = CLOSED_BOUNDARY
>>> z = mg.node_x.copy()
>>> z[[5, 13]] = z[6]  # guard nodes
>>> _ = mg.add_field('node', 'topographic__elevation', z)
>>> fr = FlowRouter(mg)
>>> cf = ChiFinder(mg)
>>> fr.run_one_step()
>>> ch_nodes = np.array([4, 8, 12, 5, 9, 13, 6, 10, 14])
>>> ch_integrand = 3.*np.ones(9, dtype=float)  # to make calc clearer
>>> chi_array = np.zeros(mg.number_of_nodes, dtype=float)
>>> cf.integrate_chi_avg_dx(ch_nodes, ch_integrand, chi_array, 0.5)
>>> chi_array.reshape(mg.shape)
array([[ 0. ,  0. ,  0. ,  0. ],
       [ 1.5,  3. ,  4.5,  0. ],
       [ 1.5,  3. ,  4.5,  0. ],
       [ 1.5,  3. ,  4.5,  0. ],
       [ 0. ,  0. ,  0. ,  0. ]])
integrate_chi_each_dx(valid_upstr_order, chi_integrand_at_nodes, chi_array)[source]

Calculates chi at each channel node by summing chi_integrand*dx.

This method accounts explicitly for spacing between each node. Method is deliberately split out for potential cythonization at a later stage. Uses a trapezium integration method.

Parameters:

valid_upstr_order : array of ints

nodes in the channel network in upstream order.

chi_integrand_at_nodes : array of floats

The value (A0/A)**concavity, in node order.

chi_array : array of floats

Array in which to store chi.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid, CLOSED_BOUNDARY
>>> from landlab.components import FlowRouter
>>> mg = RasterModelGrid((5, 4), 3.)
>>> for nodes in (mg.nodes_at_right_edge, mg.nodes_at_bottom_edge,
...               mg.nodes_at_top_edge):
...     mg.status_at_node[nodes] = CLOSED_BOUNDARY
>>> z = mg.node_x.copy()
>>> z[[5, 13]] = z[6]  # guard nodes
>>> _ = mg.add_field('node', 'topographic__elevation', z)
>>> fr = FlowRouter(mg)
>>> cf = ChiFinder(mg)
>>> fr.run_one_step()
>>> ch_nodes = np.array([4, 8, 12, 5, 9, 13, 6, 10, 14])
>>> ch_integrand = 2.*np.ones(mg.number_of_nodes,
...                           dtype=float)  # to make calc clearer
>>> chi_array = np.zeros(mg.number_of_nodes, dtype=float)
>>> cf.integrate_chi_each_dx(ch_nodes, ch_integrand, chi_array)
>>> chi_array.reshape(mg.shape)
array([[  0.        ,   0.        ,   0.        ,   0.        ],
       [  0.        ,   6.        ,  14.48528137,   0.        ],
       [  0.        ,   6.        ,  12.        ,   0.        ],
       [  0.        ,   6.        ,  14.48528137,   0.        ],
       [  0.        ,   0.        ,   0.        ,   0.        ]])
>>> from landlab.components import FastscapeEroder
>>> mg2 = RasterModelGrid((5, 5), 100.)
>>> for nodes in (mg2.nodes_at_right_edge, mg2.nodes_at_bottom_edge,
...               mg2.nodes_at_top_edge):
...     mg2.status_at_node[nodes] = CLOSED_BOUNDARY
>>> _ = mg2.add_zeros('node', 'topographic__elevation')
>>> mg2.at_node['topographic__elevation'][mg2.core_nodes] = mg2.node_x[
...     mg2.core_nodes]/1000.
>>> np.random.seed(0)
>>> mg2.at_node['topographic__elevation'][
...     mg2.core_nodes] += np.random.rand(mg2.number_of_core_nodes)
>>> fr2 = FlowRouter(mg2)
>>> sp2 = FastscapeEroder(mg2, K_sp=0.01)
>>> cf2 = ChiFinder(mg2, min_drainage_area=1., reference_concavity=0.5,
...                 use_true_dx=True)
>>> for i in range(10):
...     mg2.at_node['topographic__elevation'][mg2.core_nodes] += 10.
...     fr2.run_one_step()
...     sp2.run_one_step(1000.)
>>> fr2.run_one_step()
>>> output_array = np.zeros(25, dtype=float)
>>> cf2.integrate_chi_each_dx(mg2.at_node['flow__upstream_node_order'],
...                           np.ones(25, dtype=float),
...                           output_array)
>>> output_array.reshape(mg2.shape)
array([[   0. ,    0. ,    0.        ,    0.        ,    0. ],
       [   0. ,  100. ,  200.        ,  382.84271247,    0. ],
       [   0. ,  100. ,  241.42135624,  341.42135624,    0. ],
       [   0. ,  100. ,  200.        ,  300.        ,    0. ],
       [   0. ,    0. ,    0.        ,    0.        ,    0. ]])
masked_chi_indices

Returns a masked array version of the ‘channel__chi_index’ field. This enables easier plotting of the values with landlab.imshow_grid_at_node() or similar.

Examples

Make a topographic map with an overlay of chi values:

>>> from landlab import imshow_grid_at_node
>>> from landlab import RasterModelGrid, CLOSED_BOUNDARY
>>> from landlab.components import FlowRouter, FastscapeEroder
>>> mg = RasterModelGrid((5, 5), 100.)
>>> for nodes in (mg.nodes_at_right_edge, mg.nodes_at_bottom_edge,
...               mg.nodes_at_top_edge):
...     mg.status_at_node[nodes] = CLOSED_BOUNDARY
>>> _ = mg.add_zeros('node', 'topographic__elevation')
>>> mg.at_node['topographic__elevation'][mg.core_nodes] = mg.node_x[
...     mg.core_nodes]/1000.
>>> np.random.seed(0)
>>> mg.at_node['topographic__elevation'][
...     mg.core_nodes] += np.random.rand(mg.number_of_core_nodes)
>>> fr = FlowRouter(mg)
>>> sp = FastscapeEroder(mg, K_sp=0.01)
>>> cf = ChiFinder(mg, min_drainage_area=20000.)
>>> for i in range(10):
...     mg.at_node['topographic__elevation'][mg.core_nodes] += 10.
...     fr.run_one_step()
...     sp.run_one_step(1000.)
>>> fr.run_one_step()
>>> cf.calculate_chi()
>>> imshow_grid_at_node(mg, 'topographic__elevation',
...                     allow_colorbar=False)
>>> imshow_grid_at_node(mg, cf.masked_chi_indices,
...                     color_for_closed=None, cmap='winter')
mean_channel_node_spacing(ch_nodes)[source]

Calculates the mean spacing between all adjacent channel nodes.

Parameters:

ch_nodes : array of ints

The nodes within the defined channel network.

Returns:

mean_spacing : float (m)

The mean spacing between all nodes in the network.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid, CLOSED_BOUNDARY
>>> from landlab.components import FlowRouter
>>> mg = RasterModelGrid((5, 4), 2.)
>>> for nodes in (mg.nodes_at_right_edge, mg.nodes_at_bottom_edge,
...               mg.nodes_at_top_edge):
...     mg.status_at_node[nodes] = CLOSED_BOUNDARY
>>> z = mg.node_x.copy()
>>> z[[5, 13]] = z[6]  # guard nodes
>>> _ = mg.add_field('node', 'topographic__elevation', z)
>>> fr = FlowRouter(mg)
>>> cf = ChiFinder(mg)
>>> fr.run_one_step()
>>> ch_nodes = np.array([4, 8, 12, 5, 9, 13, 6, 10, 14])
>>> cf.mean_channel_node_spacing(ch_nodes)
2.2761423749153966
nodes_downstream_of_channel_head(channel_head)[source]

Find and return an array with nodes downstream of channel_head.

Parameters:

channel_head : int

Node ID of channel head from which to get downstream nodes.

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid, CLOSED_BOUNDARY
>>> from landlab.components import FlowRouter
>>> mg = RasterModelGrid((3, 4), 1.)
>>> for nodes in (mg.nodes_at_right_edge, mg.nodes_at_bottom_edge,
...               mg.nodes_at_top_edge):
...     mg.status_at_node[nodes] = CLOSED_BOUNDARY
>>> z = mg.add_field('node', 'topographic__elevation',
...                  mg.node_x.copy())
>>> z[4:8] = np.array([0.5, 1., 2., 0.])
>>> fr = FlowRouter(mg)
>>> fr.run_one_step()
>>> mg.at_node['flow__receiver_node']
array([ 0,  1,  2,  3,  4,  4,  5,  7,  8,  9, 10, 11])
>>> cf = ChiFinder(mg, min_drainage_area=0., reference_concavity=1.)
>>> cf.calculate_chi()
>>> cf.nodes_downstream_of_channel_head(6)
[6, 5, 4]