Gradient calculators#

Calculate gradients of quantities over links.

Gradient calculation functions#

calc_grad_at_link(node_values[, out])

Calculate gradients of node values at links.

calc_diff_at_link(node_values[, out])

Calculate differences of node values over links.

calc_aspect_at_node(grid, slope_component_tuple=None, elevs='topographic__elevation', unit='degrees', ignore_closed_nodes=True)[source]#

Get array of aspect of a surface.

Calculates at returns the aspect of a surface. Aspect is returned as radians clockwise of north, unless input parameter units is set to ‘degrees’.

If slope_component_tuple is provided, i.e., (slope_x, slope_y), the aspect will be calculated from these data.

If it is not, it will be derived from elevation data at the nodes, which can either be a string referring to a grid field (default: ‘topographic__elevation’), or an nnodes-long numpy array of the values themselves.

If ignore_closed_nodes is False, all proximal elevation values will be used in the calculation. If True, only unclosed nodes are used.

Parameters:
  • slope_component_tuple ((slope_x_array, slope_y_array) (optional)) – Tuple of components of slope in the x and y directions, defined on nodes, if already known. If not, provide elevs.

  • elevs (str or array (optional)) – Node field name or node array of elevations. If slope_component_tuple is not provided, must be set, but unused otherwise.

  • unit ({'degrees', 'radians'}) – Controls the unit that the aspect is returned as.

  • ignore_closed_nodes (bool) – If True, do not incorporate values at closed nodes into the calc.

Examples

>>> from landlab import RasterModelGrid
>>> mg = RasterModelGrid((4, 4))
>>> z = mg.node_x**2 + mg.node_y**2
>>> mg.calc_aspect_at_node(elevs=z)
array([ 225.        ,  240.16585039,  255.2796318 ,  258.69006753,
        209.83414961,  225.        ,  243.54632481,  248.77808974,
        194.7203682 ,  206.45367519,  225.        ,  231.94498651,
        191.30993247,  201.22191026,  218.05501349,  225.        ])
>>> z = z.max() - z
>>> mg.calc_aspect_at_node(elevs=z)
array([ 45.        ,  60.16585039,  75.2796318 ,  78.69006753,
        29.83414961,  45.        ,  63.54632481,  68.77808974,
        14.7203682 ,  26.45367519,  45.        ,  51.94498651,
        11.30993247,  21.22191026,  38.05501349,  45.        ])
>>> mg = RasterModelGrid((4, 4), xy_spacing=(3.0, 2.0))
>>> z = mg.node_x**2 + mg.node_y**2
>>> mg.calc_aspect_at_node(elevs=z)
array([ 236.30993247,  247.52001262,  259.97326008,  262.40535663,
        220.75264634,  234.41577266,  251.13402374,  255.29210302,
        201.54258265,  215.47930877,  235.73541937,  242.24162456,
        196.69924423,  209.43534223,  229.19345757,  236.30993247])

Note that a small amount of asymmetry arises at the grid edges due to the “missing” nodes beyond the edge of the grid.

Calculate differences of node values over links.

Calculates the difference in quantity node_values at each link in the grid.

Parameters:
  • node_values (ndarray or field name) – Values at grid nodes.

  • out (ndarray, optional) – Buffer to hold the result.

Returns:

Differences across links.

Return type:

ndarray

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> rmg = RasterModelGrid((3, 3))
>>> z = np.zeros(9)
>>> z[4] = 1.0
>>> rmg.calc_diff_at_link(z)
array([ 0.,  0.,  0.,  1.,  0.,  1., -1.,  0., -1.,  0.,  0.,  0.])

Calculate gradients of node values at links.

Calculates the gradient in node_values at each link in the grid, returning an array of length number_of_links.

Parameters:
  • node_values (ndarray or field name (x number of nodes)) – Values at grid nodes.

  • out (ndarray, optional (x number of links)) – Buffer to hold the result.

Returns:

Gradients across active links.

Return type:

ndarray

Examples

>>> from landlab import RasterModelGrid
>>> rg = RasterModelGrid((3, 4), xy_spacing=10.0)
>>> z = rg.add_zeros("topographic__elevation", at="node")
>>> z[5] = 50.0
>>> z[6] = 36.0
>>> calc_grad_at_link(rg, z)  # there are 17 links
array([ 0. ,  0. ,  0. ,  0. ,  5. ,  3.6,  0. ,  5. , -1.4, -3.6,  0. ,
       -5. , -3.6,  0. ,  0. ,  0. ,  0. ])
>>> from landlab import HexModelGrid
>>> hg = HexModelGrid((3, 3), spacing=10.0)
>>> z = hg.add_zeros("topographic__elevation", at="node", clobber=True)
>>> z[4] = 50.0
>>> z[5] = 36.0
>>> calc_grad_at_link(hg, z)  # there are 11 faces
array([ 0. ,  0. ,  0. ,  5. ,  5. ,  3.6,  3.6,  0. ,  5. , -1.4, -3.6,
        0. , -5. , -5. , -3.6, -3.6,  0. ,  0. ,  0. ])
calc_grad_at_patch(grid, elevs='topographic__elevation', ignore_closed_nodes=True, unit_normal=None, slope_magnitude=None)[source]#

Calculate the components of the gradient at each patch.

If ignore_closed_nodes is True, closed nodes do not affect gradient calculations. If a closed node is present in a patch, the patch gradient is set to zero in both x and y directions.

Parameters:
  • elevs (str or ndarray, optional) – Field name or array of node values.

  • ignore_closed_nodes (bool) – If True, do not incorporate values at closed nodes into the calc.

  • unit_normal (array with shape (num_patches, 3) (optional)) – The unit normal vector to each patch, if already known.

  • slope_magnitude (array with size num_patches (optional)) – The slope of each patch, if already known.

Returns:

gradient_tuple – Len-2 tuple of arrays giving components of gradient in the x and y directions, in the units of units.

Return type:

(x_component_at_patch, y_component_at_patch)

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> mg = RasterModelGrid((4, 5))
>>> z = mg.node_y
>>> (x_grad, y_grad) = mg.calc_grad_at_patch(elevs=z)
>>> np.allclose(y_grad, np.pi / 4.0)
True
>>> np.allclose(x_grad, 0.0)
True
calc_slope_at_node(grid, elevs='topographic__elevation', method='patch_mean', ignore_closed_nodes=True, return_components=False, **kwds)[source]#

Array of slopes at nodes, averaged over neighboring patches.

Produces a value for node slope (i.e., mean gradient magnitude) at each node in a manner analogous to a GIS-style slope map. It averages the gradient on each of the patches surrounding the node, creating a value for node slope that better incorporates nonlocal elevation information. Directional information can still be returned through use of the return_components keyword.

Note that under these definitions, it is not always true that:

mag, cmp = mg.calc_slope_at_node(z)
mag**2 == cmp[0] ** 2 + cmp[1] ** 2  # not always true

If ignore_closed_nodes is False, all proximal elevation values will be used in the calculation. If True, only unclosed nodes are used.

Parameters:
  • elevs (str or ndarray, optional) – Field name or array of node values.

  • method ({'patch_mean', 'Horn'}) – By equivalence to the raster version, 'patch_mean' returns a scalar mean on the patches; 'Horn' returns a vector mean on the patches.

  • ignore_closed_nodes (bool) – If True, do not incorporate values at closed nodes into the calculation.

  • return_components (bool) – If True, return a tuple, (array_of_magnitude, (array_of_slope_x_radians, array_of_slope_y_radians)). If False, return an array of floats of the slope magnitude.

Returns:

If return_components, returns (array_of_magnitude, (array_of_slope_x_radians, array_of_slope_y_radians)). If not return_components, returns an array of slope magnitudes.

Return type:

float array or length-2 tuple of float arrays

Examples

>>> import numpy as np
>>> from landlab import RadialModelGrid, RasterModelGrid
>>> mg = RasterModelGrid((4, 5))
>>> z = mg.node_x
>>> slopes = mg.calc_slope_at_node(elevs=z)
>>> np.allclose(slopes, 45.0 / 180.0 * np.pi)
True
>>> mg = RasterModelGrid((4, 5))
>>> z = -mg.node_y
>>> slope_mag, cmp = mg.calc_slope_at_node(elevs=z, return_components=True)
>>> np.allclose(slope_mag, np.pi / 4.0)
True
>>> np.allclose(cmp[0], 0.0)
True
>>> np.allclose(cmp[1], -np.pi / 4.0)
True
>>> mg = RadialModelGrid(n_rings=3)
>>> z = mg.radius_at_node
>>> slope_at_node = np.round(mg.calc_slope_at_node(elevs=z), decimals=5)
>>> nodes_at_ring = [
...     np.where(np.isclose(mg.radius_at_node, radius)) for radius in range(3)
... ]
>>> slope_at_node[nodes_at_ring[0]]
array([ 0.85707])
>>> slope_at_node[nodes_at_ring[1]]
array([ 0.79417,  0.79417,  0.79417,  0.79417,  0.79417,  0.79417])
>>> slope_at_node[nodes_at_ring[2]]
array([ 0.77542,  0.78453,  0.78453,  0.77542,  0.77542,  0.78453,
        0.78453,  0.77542,  0.77542,  0.78453,  0.78453,  0.77542])
calc_slope_at_patch(grid, elevs='topographic__elevation', ignore_closed_nodes=True, unit_normal=None)[source]#

Calculate the slope (positive magnitude of gradient) at patches.

If ignore_closed_nodes is True, closed nodes do not affect slope calculations. If a closed node is present in a patch, the patch slope is set to zero.

Parameters:
  • elevs (str or ndarray, optional) – Field name or array of node values.

  • ignore_closed_nodes (bool) – If True, do not incorporate values at closed nodes into the calc.

  • unit_normal (array with shape (num_patches, 3) (optional)) – The unit normal vector to each patch, if already known.

Returns:

slopes_at_patch – The slope (positive gradient magnitude) of each patch.

Return type:

n_patches-long array

Examples

>>> import numpy as np
>>> from landlab import RasterModelGrid
>>> mg = RasterModelGrid((4, 5))
>>> z = mg.node_x
>>> S = mg.calc_slope_at_patch(elevs=z)
>>> S.size == mg.number_of_patches
True
>>> np.allclose(S, np.pi / 4.0)
True
calc_unit_normal_at_patch(grid, elevs='topographic__elevation')[source]#

Calculate and return the unit normal vector <a, b, c> to a patch.

Parameters:

elevs (str or ndarray, optional) – Field name or array of node values.

Returns:

nhat – The unit normal vector <a, b, c> to each patch.

Return type:

num-patches x length-3 array

Examples

>>> from landlab import HexModelGrid
>>> mg = HexModelGrid((3, 3))
>>> z = mg.node_x * 3.0 / 4.0
>>> mg.calc_unit_normal_at_patch(z)
array([[-0.6,  0. ,  0.8],
       [-0.6,  0. ,  0.8],
       [-0.6,  0. ,  0.8],
       [-0.6,  0. ,  0.8],
       [-0.6,  0. ,  0.8],
       [-0.6,  0. ,  0.8],
       [-0.6,  0. ,  0.8],
       [-0.6,  0. ,  0.8],
       [-0.6,  0. ,  0.8],
       [-0.6,  0. ,  0.8]])