"""Function to calculate node statistics in a moving window."""
import numpy as np
from landlab import FieldError
[docs]
def calculate_window_statistic(
grid, field, func, search_radius, calc_on_closed_nodes=True, **kwargs
):
"""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 :class:`~.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 :func:`~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 : ndarray
Output array containing the calculated values of the statistic. Same
length as input field.
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])
"""
if field not in grid.at_node:
raise FieldError(f"A {field} field is required at the nodes of the input grid.")
# Create output array
output = np.zeros(grid.number_of_nodes)
# Create arrays of x and y coords for input to "distance to point' calc
x_coord = grid.x_of_node
y_coord = grid.y_of_node
nodes_in_loop = grid.nodes.flatten()
nodes_to_include = np.ones(grid.number_of_nodes, dtype=bool)
if calc_on_closed_nodes is False:
closed_nodes = grid.status_at_node == grid.BC_NODE_IS_CLOSED
nodes_in_loop = nodes_in_loop[~closed_nodes]
nodes_to_include[closed_nodes] = False
output[closed_nodes] = np.nan
# Calculate "dist to point" then local value at nodes within window.
for node in nodes_in_loop:
node_dist_to_point = grid.calc_distances_of_nodes_to_point(
(x_coord[node], y_coord[node])
)
nodes_in_window = np.all(
[node_dist_to_point <= search_radius, nodes_to_include], 0
)
values_in_window = grid.at_node[field][nodes_in_window]
output[node] = func(values_in_window, **kwargs)
return output