# Source code for landlab.utils.structured_grid

```#! /usr/bin/env python
"""Utility functions for structured grid of elements with four neighbors."""

import contextlib
import itertools

import numpy as np

from ..core.utils import as_id_array
from ..grid.nodestatus import NodeStatus

[docs]
def node_count(shape):
"""Total number of nodes.

The total number of nodes in a structured grid with dimensions given
by the tuple, *shape*. Where *shape* is the number of node rows and
node columns.

>>> from landlab.utils.structured_grid import node_count
>>> node_count((3, 4))
12
"""
assert len(shape) == 2
return shape * shape

[docs]
def interior_node_count(shape):
"""Number of interior nodes.

Return the count of the number of interior nodes of a structured grid
of dimensions, *shape*.

>>> from landlab.utils.structured_grid import node_count, interior_node_count
>>> node_count((2, 4))
8
>>> interior_node_count((2, 4))
0
>>> interior_node_count((1, 4))
0
>>> interior_node_count((3, 4))
2
"""
assert len(shape) == 2
if np.min(shape) > 2:
return (shape - 2) * (shape - 2)
else:
return 0

[docs]
def cell_count(shape):
"""Total number of cells.

The total number of cells in a structured grid with dimensions, *shape*.
Where *shape* is a tuple that gives the dimensions of the grid as number
of rows of nodes followed by number of columns of nodes.

>>> from landlab.utils.structured_grid import cell_count
>>> cell_count((3, 4))
2
>>> cell_count((1, 4))
0
"""
assert len(shape) == 2
if np.min(shape) > 2:
return (shape - 2) * (shape - 2)
else:
return 0

[docs]
def active_cell_count(shape):
"""Number of active cells.

Number of active cells. By default, all cells are active so this is
the same as cell_count. (active = core+open boundary)
"""
return cell_count(shape)

[docs]
def core_cell_count(shape):
"""Number of core cells.

Number of core cells. By default, all cells are core so this is the
same as cell_count.
"""
return cell_count(shape)

[docs]

Number of active links in a structured grid with dimensions, *shape*.
A link is active if it connects to at least one active node.

7
0
7
"""
assert len(shape) == 2
if np.min(shape) > 2:
return 2 * shape * shape - 3 * (shape + shape) + 4
else:
return 0

[docs]

Total (active and inactive) number of links in a structured grid with
the total of the two.

7
"""
assert len(shape) == 2
return shape * (shape - 1) + shape * (shape - 1)

[docs]
assert len(shape) == 2
return (shape - 1) * shape

[docs]
assert len(shape) == 2
return shape * (shape - 1)

[docs]
def perimeter_node_count(shape):
"""Number of perimeter nodes.

Number of nodes that are on the perimeter of a structured grid with
dimensions, *shape*, and thus boundary nodes.

Examples
--------
>>> from landlab.utils.structured_grid import perimeter_node_count
>>> perimeter_node_count((3, 4))
10
"""
assert len(shape) == 2
return 2 * (shape - 2) + 2 * (shape - 2) + 4

[docs]
def interior_cell_count(shape):
"""Number of interior cells.

Number of interior cells. Since cells are only defined on interior
nodes, this is the same as cell_count.
"""
return cell_count(shape)

[docs]
def face_count(shape):
"""Total number of faces.

Total number of faces in a structured grid with dimensions, *shape*. Each
cell has four faces, and shared faces only count once.

Examples
--------
>>> from landlab.utils.structured_grid import face_count
>>> face_count((3, 4))
7
"""
assert len(shape) == 2
if np.min(shape) > 2:
return (shape - 1) * (shape - 2) + (shape - 2) * (shape - 1)
else:
return 0

[docs]
def active_face_count(shape):
"""Number of active faces.

Total number of active faces in a structured grid with dimensions,
*shape*. Each cell has four faces, and shared faces only count once.
An active face is one that has a corresponing active link.

>>> from landlab.utils.structured_grid import active_face_count
>>> active_face_count((3, 4))
7
"""
return face_count(shape)

[docs]
def top_index_iter(shape):
"""Iterator for the top boundary indices of a structured grid."""
return range(shape * (shape - 1), shape * shape)

[docs]
def bottom_index_iter(shape):
"""Iterator for the bottom boundary indices of a structured grid."""
return range(0, shape)

[docs]
def left_index_iter(shape):
"""Iterator for the left boundary indices of a structured grid."""
return range(0, shape * shape, shape)

[docs]
def right_index_iter(shape):
"""Iterator for the right boundary indices of a structured grid."""
return range(shape - 1, shape * shape, shape)

[docs]
def left_right_iter(shape, *args):
"""Iterator for the left and right boundary indices of a structured grid.

This iterates over the indices in order rather than iterating all of
the left boundary and then all of the right boundary.

Examples
--------
>>> import numpy as np
>>> from landlab.utils.structured_grid import left_right_iter
>>> np.fromiter(left_right_iter((4, 3)), dtype=int)
array([ 0,  2,  3,  5,  6,  8,  9, 11])
>>> np.fromiter(left_right_iter((4, 3), 2), dtype=int)
array([0, 2, 3, 5])
>>> np.fromiter(left_right_iter((4, 3), 2, 4), dtype=int)
array([ 6,  8,  9, 11])
>>> np.fromiter(left_right_iter((4, 3), 1, 4, 2), dtype=int)
array([ 3,  5,  9, 11])
"""
if len(args) == 0:
iter_rows = range(0, shape, 1)
elif len(args) == 1:
iter_rows = range(0, args, 1)
elif len(args) == 2:
iter_rows = range(args, args, 1)
elif len(args) == 3:
iter_rows = range(args, args, args)

for row in iter_rows:
yield row * shape
yield row * shape + shape - 1

[docs]
def bottom_top_iter(shape):
"""Iterator for bottom then top indices of a structured grid.

Examples
--------
>>> import numpy as np
>>> from landlab.utils.structured_grid import bottom_top_iter
>>> np.fromiter(bottom_top_iter((4, 3)), dtype=int)
array([ 0,  1,  2,  9, 10, 11])
"""
return itertools.chain(bottom_index_iter(shape), top_index_iter(shape))

[docs]
def perimeter_iter(shape):
"""Iterator for perimeter nodes.

Iterates over all of the perimeter node indices of a structured grid in
order.

Examples
--------
>>> import numpy as np
>>> from landlab.utils.structured_grid import perimeter_iter
>>> np.fromiter(perimeter_iter((4, 3)), dtype=int)
array([ 0,  1,  2,  3,  5,  6,  8,  9, 10, 11])
"""
return itertools.chain(
bottom_index_iter(shape),
left_right_iter(shape, 1, shape - 1),
top_index_iter(shape),
)

[docs]
def perimeter_nodes(shape):
"""Array of perimeter nodes.

An array of the indices of the perimeter nodes of a structured grid.

Examples
--------
>>> from landlab.utils.structured_grid import perimeter_nodes
>>> perimeter_nodes((3, 4))
array([ 0,  1,  2,  3,  4,  7,  8,  9, 10, 11])
"""
return np.fromiter(perimeter_iter(shape), dtype=int)

[docs]
def corners(shape):
"""Array of the indices of the grid corner nodes."""
return np.array(
[0, shape - 1, shape * (shape - 1), shape * shape - 1]
)

[docs]
def bottom_edge_node_ids(shape):
"""Array of nodes on the bottom edge."""
return np.fromiter(bottom_index_iter(shape), dtype=int)

[docs]
def top_edge_node_ids(shape):
"""Array of nodes on the top edge."""
return np.fromiter(top_index_iter(shape), dtype=int)

[docs]
def left_edge_node_ids(shape):
"""Array of nodes on the left edge."""
return np.fromiter(left_index_iter(shape), dtype=int)

[docs]
def right_edge_node_ids(shape):
"""Array of nodes on the right edge."""
return np.fromiter(right_index_iter(shape), dtype=int)

[docs]
def interior_iter(shape):
"""Iterator for the interior nodes of a structured grid.

Examples
--------
>>> import numpy as np
>>> from landlab.utils.structured_grid import interior_iter
>>> np.fromiter(interior_iter((4, 3)), dtype=int)
array([4, 7])
"""
interiors = []
interiors_per_row = shape - 2
for row in range(shape + 1, shape * (shape - 1), shape):
interiors.append(range(row, row + interiors_per_row))
return itertools.chain(*interiors)

[docs]
def interior_nodes(shape):
"""Array of interior nodes."""
return np.fromiter(interior_iter(shape), dtype=int)

[docs]
def node_coords(shape, *args):
"""Get node x and y coordinates.

Get x, y coordinates for nodes in a structured grid with dimensions,
*shape*. Use the optional argument *spacing* to give the spacing in each
dimension, and *origin* the start of the coordinates in each dimension.

Parameters
----------
shape : tuple of int
Number of node rows and columns.
spacing : tuple, optional
Row and column spacing.
origin : tuple, optional
Coordinate of lower-left node.

Examples
--------
>>> from landlab.utils.structured_grid import node_coords
>>> (cols, rows) = node_coords((3, 2))
>>> rows
array([ 0.,  0.,  1.,  1.,  2.,  2.])
>>> cols
array([ 0.,  1.,  0.,  1.,  0.,  1.])
"""
try:
spacing = args
except IndexError:
spacing = np.ones(len(shape), dtype=float)
else:
assert len(spacing) == len(shape)

try:
origin = args
except IndexError:
origin = np.zeros(len(shape), dtype=float)
else:
assert len(origin) == len(origin)

node_count_ = np.prod(shape)

row_y = np.arange(shape) * spacing + origin
col_x = np.arange(shape) * spacing + origin

(node_x, node_y) = np.meshgrid(col_x, row_y)

node_x.shape = (node_count_,)
node_y.shape = (node_count_,)

return (node_x, node_y)

[docs]
def node_at_cell(shape):
"""Array of nodes at cells.

Indices of the nodes belonging to each cell.

Examples
--------
>>> from landlab.utils.structured_grid import node_at_cell
>>> node_at_cell((4, 3))
array([4, 7])
"""
node_ids = np.arange(node_count(shape))
node_ids.shape = shape

cell_node = node_ids[1:-1, 1:-1].copy()
cell_node.shape = ((shape - 2) * (shape - 2),)

return cell_node

[docs]
"""Array of nodes at each end of links."""
node_ids = np.arange(np.prod(shape))
node_ids.shape = shape

[docs]

[docs]

[docs]
"""Array of nodes at the end of links."""

[docs]
"""Array of nodes at the start of links."""

[docs]
"""Array of faces associated with links.

Returns an array that maps link ids to face ids. For inactive links,
which do not have associated faces, set their ids to
*inactive_link_index*. Use the *actives* keyword to specify an array that
contains the ids of all active links in the grid. The default assumes
that only the perimeter nodes are inactive.

Examples
--------
>>> faces
array([-1,  0,  1, -1, -1,  2,  3,
-1, -1, -1, -1,  4,  5,  6, -1, -1, -1])
"""
if actives is None:

[docs]
def status_at_node(shape, boundary_status=NodeStatus.FIXED_VALUE):
"""Array of the statuses of nodes.

The statuses of the nodes in a structured grid with dimensions,
*shape*. Use the *boundary_status* keyword to specify the status of
the top, bottom, left and right boundary nodes.
"""
status = np.empty(np.prod(shape), dtype=np.int8)

status[interior_nodes(shape)] = NodeStatus.CORE
status[perimeter_nodes(shape)] = boundary_status

return status

[docs]

Return the link IDs for links that are *active* in a structured grid of
quadrilaterals. Use the *node_status_array* keyword to specify the status
for each of the grid's nodes. If not given, each of the perimeter nodes is
assumed to be `NodeStatus.FIXED_VALUE`.

Use the *link_nodes* keyword to provide, as a tuple of arrays, that give
the *from-node* and the *to-node* for each for each link in the grid.

Parameters
----------
shape : tuple
Shape of grid as number of node rows and columns.
node_status_array : array_like, optional
Status of each grid node.

Examples
--------
Because, by default, the perimeter nodes are `NodeStatus.FIXED_VALUE` nodes,
only links attached to the interior nodes are *active*.

>>> from landlab.grid.nodestatus import NodeStatus
array([ 1,  2,  5,  6, 11, 12, 13])

If all the perimeter nodes `NodeStatus.CLOSED` nodes, the only active link
is between the two core nodes.

>>> node_status = np.ones(3 * 4) * NodeStatus.CLOSED
>>> node_status[5:7] = NodeStatus.CORE
array()

You can also provide a list of all the *from_nodes* and *to_nodes* for
the grid. The following describes a grid with only a single link (between
nodes 5 and 6).

array()
"""
if node_status_array is None:
node_status_array = status_at_node(shape)

else:

(from_node_status == NodeStatus.CORE) & ~(to_node_status == NodeStatus.CLOSED)
) | ((to_node_status == NodeStatus.CORE) & ~(from_node_status == NodeStatus.CLOSED))

[docs]
def active_face_index(shape):
"""Array of face indices."""
return np.arange(active_face_count(shape))

[docs]

[docs]

[docs]
"""Array of active links entering nodes."""
)

[docs]
"""Array of active links entering nodes.

Finds and returns the link IDs of active links coming in to each node

Parameters
----------
shape : 2-element tuple of ints
Number of rows and columns in the grid
node_status (optional) : numpy array of bool (x # of nodes)
False where node is a closed boundary; True elsewhere

Returns
-------
2d numpy array of int (2 x number of grid nodes)

Examples
--------
array([[-1, -1, -1, -1, -1,  4,  5, -1, -1, 11, 12, -1],
[-1, -1, -1, -1, -1,  7,  8,  9, -1, -1, -1, -1]])

Notes
-----
There are at most two inlinks for each node. The first row in the returned
array gives the ID of the vertical incoming link from below (south), or -1
if there is none. The second row gives the link ID of the horizontal link
coming in from the left (or -1).
"""
(
)
)

[docs]
"""Array of active links leaving nodes."""
)

[docs]
"""Array of active links leaving nodes.

Finds and returns the link IDs of active links going out of each node
(that is, active links for which the node is the link tail).

Parameters
----------
shape : 2-element tuple of ints
Number of rows and columns in the grid
node_status (optional) : numpy array of bool (x # of nodes)
False where node is a closed boundary; True elsewhere

Returns
-------
2d numpy array of int (2 x number of grid nodes)

Examples
--------
array([[-1,  4,  5, -1, -1, 11, 12, -1, -1, -1, -1, -1],
[-1, -1, -1, -1,  7,  8,  9, -1, -1, -1, -1, -1]])

Notes
-----
There are at most two outlinks for each node. The first row in the returned
array gives the ID of the vertical outgoing link to above (north), or -1
if there is none. The second row gives the link ID of the horizontal link
going out to the right (east) (or -1).
"""
(
)
)

[docs]
link_ids = np.empty((shape - 1, shape), dtype=int)
num_links_per_row = (2 * shape) - 1
for r in range(shape - 1):
link_ids[r, :] = (shape - 1) + (r * num_links_per_row) + np.arange(shape)

[docs]
link_ids = np.empty((shape, shape - 1), dtype=int)
num_links_per_row = (2 * shape) - 1
for r in range(shape):

[docs]
"""Number of active links oriented vertically."""
if node_status is not None:
is_inactive = node_status == 0
is_inactive.shape = shape

return (shape - 1) * (shape - 2) - np.sum(inactive_links.flat)
else:
return (shape - 1) * (shape - 2)

[docs]
"""Number of active links oriented horizontally."""
if node_status is not None:
is_inactive = node_status == 0
is_inactive.shape = shape

return (shape - 2) * (shape - 1) - np.sum(inactive_links.flat)
else:
return (shape - 2) * (shape - 1)

[docs]

Creates and returns a boolean 2D array dimensioned as the number of
vertical links in the grid, not including the left and right boundaries.

Parameters
----------
shape : 2-element tuple of ints
Number of rows and columns in the grid
node_status : numpy array of bool (x # of nodes)
False where node is a closed boundary; True elsewhere

Returns
-------
(NR-1,NC-2) array of bool (NR=# of rows, NC=# of columns)
Flags indicating whether the corresponding vertical link is inactive

Examples
--------
>>> import numpy as np
>>> ns = np.ones(12, dtype=bool)  # case of no closed boundary nodes
array([[False, False],
[False, False]], dtype=bool)
>>> ns = False  # node 2 is a closed boundary
array([[False,  True],
[False, False]], dtype=bool)
>>> ns = False  # node 9 is also a closed boundary
array([[False,  True],
[ True, False]], dtype=bool)
"""
# Create a 2D boolean matrix indicating whether NODES are closed boundaries
# GT thinks this should be False, not 0
is_closed_node = node_status == 0
is_closed_node.shape = shape

inactive_outlinks = is_closed_node[:-1, 1:-1]  # middle cols, all but top row
# middle cols, all but bottom row
# if either node is closed, the link is inactive

[docs]
is_inactive = node_status == 0
is_inactive.shape = shape

#    link_ids.shape = (shape - 1, shape - 2)

[docs]
"""Array of active links oriented vertically."""
if node_status is None:
# link_ids.shape = (shape - 1, shape - 2)
else:

link_ids.shape = (shape - 1, shape - 2)

[docs]
"""Array of active links oriented vertically.

Returns the link IDs of vertical active links as an (R-1) x (C-2) array.

Parameters
----------
shape : 2-element tuple of int
number of rows and columns in grid
node_status (optional) : 1d numpy array (x number of nodes) of bool
False where node is a closed boundary, True otherwise

Returns
-------
2d numpy array of int
left and right grid edges. If a vertical link is inactive, its ID is
given as -1.

Examples
--------
array([[ 4,  5],
[11, 12]])
>>> ns = np.ones(12, dtype=bool)
>>> ns = False
>>> ns = False
array([[-1,  5],
[11, -1]])

Notes
-----
"""
link_ids = np.empty((shape - 1, shape - 2), dtype=int)
num_links_per_row = (2 * shape) - 1
for r in range(shape - 1):
link_ids[r, :] = shape + (r * num_links_per_row) + np.arange(shape - 2)

if node_status is not None:

[docs]
"""Array of active links oriented horizontally."""
if node_status is None:
dtype=int,
)
else:

)
link_ids.shape = (shape - 2, shape - 1)

[docs]
"""Array of active links oriented horizontally.

Returns the link IDs of horizontal active links as an (R-2) x (C-1) array.

Parameters
----------
shape : 2-element tuple of int
number of rows and columns in grid
node_status (optional) : 1d numpy array (x number of nodes) of bool
False where node is a closed boundary, True otherwise

Returns
-------
2d numpy array of int
top and bottom grid edges. If a horizontal link is inactive, its ID is
given as -1.

Examples
--------
array([[7, 8, 9]])
>>> ns = np.ones(12, dtype=bool)
>>> ns = False
>>> ns = False
array([[-1,  8, -1]])

Notes
-----
"""
link_ids = np.empty((shape - 2, shape - 1), dtype=int)
num_links_per_row = (2 * shape) - 1
for r in range(shape - 2):
link_ids[r, :] = ((r + 1) * num_links_per_row) + np.arange(shape - 1)
if node_status is not None:

[docs]
"""Array of links pointing to the west.

Examples
--------
array([[-1,  0,  1,  2],
[-1,  7,  8,  9],
[-1, 14, 15, 16]])
"""
link_ids.shape = (shape, shape - 1)

[docs]
"""Array of links pointing to the north.

Examples
--------
array([[ 3,  4,  5,  6],
[10, 11, 12, 13],
[-1, -1, -1, -1]])
"""
link_ids.shape = (shape - 1, shape)

[docs]
"""Array of links pointing to the the south.

Examples
--------
array([[-1, -1, -1, -1],
[ 3,  4,  5,  6],
[10, 11, 12, 13]])
"""
link_ids.shape = (shape - 1, shape)

[docs]
"""Array of links pointing to the the east.

Examples
--------
array([[ 0,  1,  2, -1],
[ 7,  8,  9, -1],
[14, 15, 16, -1]])
"""
link_ids.shape = (shape, shape - 1)

[docs]
"""Array of active links pointing to the the north."""
with contextlib.suppress(ValueError):
links.shape = (shape - 1, shape - 2)

[docs]
"""Array of active links pointing to the the north.

array([[-1,  4,  5, -1],
[-1, 11, 12, -1],
[-1, -1, -1, -1]])
"""
with contextlib.suppress(ValueError):
links.shape = (shape - 1, shape - 2)

[docs]
"""Array of active links pointing to the the south."""
links.shape = (shape - 1, shape - 2)

[docs]
"""Array of active links pointing to the the south.

Finds and returns link IDs of active links that enter each node from the
south (bottom), or -1 where no such active link exists.

Parameters
----------
shape : 2-element tuple of int
number of rows and columns in grid
node_status (optional) : 1d numpy array of bool
False where node is a closed boundary, True otherwise

Returns
-------
2d numpy array of int
Link ID of active link connecting to a node from the south, or -1

Examples
--------
array([[-1, -1, -1, -1],
[-1,  4,  5, -1],
[-1, 11, 12, -1]])

Notes
-----
"""

[docs]
"""Array of active links pointing to the the west.

Examples
--------
array([[-1, -1, -1, -1],
[-1,  4,  5,  6],
[-1, -1, -1, -1]])
"""
with contextlib.suppress(ValueError):
shape, node_status=node_status
)

[docs]
"""Array of active links pointing to the the west.

Examples
--------
array([[-1, -1, -1, -1],
[-1,  7,  8,  9],
[-1, -1, -1, -1]])
"""
shape, node_status=node_status
)

[docs]
"""Array of active links pointing to the the east.

Examples
--------
array([[-1, -1, -1, -1],
[ 4,  5,  6, -1],
[-1, -1, -1, -1]])
"""
with contextlib.suppress(ValueError):
shape, node_status=node_status
)

[docs]
"""Array of active links pointing to the the east.

Examples
--------
array([[-1, -1, -1, -1],
[ 7,  8,  9, -1],
[-1, -1, -1, -1]])
"""
shape, node_status=node_status
)

[docs]
"""Number of links leaving each node."""

[docs]
"""Number of links entering each node."""

[docs]
"""Number of active links leaving each node."""

[docs]
"""Number of active links entering each node."""

[docs]
"""Create a matrix of links leaving each node."""
if return_count:
else:

[docs]
"""Create a matrix of links entering each node."""
if return_count:
else:

[docs]
"""Create a matrix of active links leaving each node."""
if return_count:
else:

[docs]
"""Create a matrix of active links leaving each node.

Return the link IDs of the active links that leave each node of a grid. The
shape of the returned array is (2, *N*) where *N* is the number of nodes in
the grid. The first row contains the link ID exiting the node to the
top, and the second row the link exiting the node to the right.

Use the *return_count* keyword to, in addition to the link IDs, return the
number of active links attached to each grid node.

Use the *node_status_array* keyword to specify the status for each of the
grid's nodes. If not given, each of the perimeter nodes is assumed to be
`NodeStatus.FIXED_VALUE`.

Parameters
----------
shape : tuple
Shape of the structured grid
node_status : array_like, optional
Status of each node in the grid.
return_count : boolean, optional
If `True`, also return an array of active link counts per node.

Returns
-------
Active link IDs for each node.
count : ndarray
Number of active links per node.

Examples
--------
Get the active link IDs for a grid of 3 nodes by 4 nodes. The first row
lists links entering nodes from the bottom, and the second links entering
from the left.

array([[-1,  4,  5, -1, -1, 11, 12, -1, -1, -1, -1, -1],
[-1, -1, -1, -1,  7,  8,  9, -1, -1, -1, -1, -1]])
>>> _, count = setup_active_outlink_matrix2((3, 4))
>>> count
array([0, 1, 1, 0, 1, 2, 2, 0, 0, 0, 0, 0])
"""
if return_count:
else:

[docs]
"""Create a matrix of active links entering each node.

Return the IDs of the active links that enter each node of a grid. The
shape of the returned array is (2, *N*) where *N* is the number of nodes
in the grid. The first row contains the link ID entering the node from the
bottom, and the second row the link entering the node from the left.

Use the *return_count* keyword to, in addition to the link IDs, return
the number of active links attached to each grid node.

Use the *node_status_array* keyword to specify the status for each of
the grid's nodes. If not given, each of the perimeter nodes is assumed
to be `NodeStatus.FIXED_VALUE`.

Parameters
----------
shape : tuple
Shape of the structured grid
node_status : array_like, optional
Status of each node in the grid.
return_count : boolean, optional
If `True`, also return an array of active link counts per node.

Returns
-------
Active link IDs for each node.
count : ndarray
Number of active links per node.

Examples
--------
Get the active link IDs for a grid of 3 nodes by 4 nodes. The first row
list links entering nodes from the bottom, and the second links entering
from the left.

array([[-1, -1, -1, -1, -1,  0,  1, -1, -1,  2,  3, -1],
[-1, -1, -1, -1, -1,  4,  5,  6, -1, -1, -1, -1]])
>>> _, count = setup_active_inlink_matrix((3, 4))
>>> count
array([0, 0, 0, 0, 0, 2, 2, 1, 0, 1, 1, 0])
"""
if return_count:
else:

[docs]
"""Create a matrix of active links entering each node.

Return the link IDs of the active links that enter each node of a grid. The
shape of the returned array is (2, *N*) where *N* is the number of nodes in
the grid. The first row contains the link ID entering the node from the
bottom, and the second row the link entering the node from the left.

Use the *return_count* keyword to, in addition to the link IDs, return the
number of active links attached to each grid node.

Use the *node_status_array* keyword to specify the status for each of the
grid's nodes. If not given, each of the perimeter nodes is assumed to be
`NodeStatus.FIXED_VALUE`.

Parameters
----------
shape : tuple
Shape of the structured grid
node_status : array_like, optional
Status of each node in the grid.
return_count : boolean, optional
If `True`, also return an array of active link counts per node.

Returns
-------
Active link IDs for each node.
count : ndarray
Number of active links per node.

Examples
--------
Get the active link IDs for a grid of 3 nodes by 4 nodes. The first row
lists links entering nodes from the bottom, and the second links entering
from the left.

array([[-1, -1, -1, -1, -1,  4,  5, -1, -1, 11, 12, -1],
[-1, -1, -1, -1, -1,  7,  8,  9, -1, -1, -1, -1]])
>>> _, count = setup_active_inlink_matrix2((3, 4))
>>> count
array([0, 0, 0, 0, 0, 2, 2, 1, 0, 1, 1, 0])
"""
if return_count:
else:

[docs]
"""Array of links with a halo of no-data values.

Examples
--------
>>> from landlab.utils.structured_grid import node_index_with_halo
>>> node_index_with_halo((2, 3), halo_indices=-1)
array([[-1, -1, -1, -1, -1],
[-1,  0,  1,  2, -1],
[-1,  3,  4,  5, -1],
[-1, -1, -1, -1, -1]])
"""
shape_with_halo = np.array(shape) + 2

ids = np.empty(shape_with_halo, dtype=int)

(interiors, boundaries) = (
interior_nodes(shape_with_halo),
perimeter_nodes(shape_with_halo),
)

ids.flat[interiors] = range(interior_node_count(shape_with_halo))
ids.flat[boundaries] = halo_indices

return ids

[docs]
"""Array of cells with a halo of no-data values.

Examples
--------
>>> from landlab.utils.structured_grid import cell_index_with_halo
>>> cell_index_with_halo((2, 3), halo_indices=-1)
array([[-1, -1, -1, -1, -1],
[-1,  0,  1,  2, -1],
[-1,  3,  4,  5, -1],
[-1, -1, -1, -1, -1]])

>>> cell_index_with_halo((2, 3), halo_indices=-1, inactive_indices=-1)
array([[-1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1]])
"""
ids = node_index_with_halo(shape, halo_indices=halo_indices)
if inactive_indices is not None:
ids[:, (1, -2)] = inactive_indices
ids[(1, -2), :] = inactive_indices

return ids

def _neighbor_node_ids(ids_with_halo):
"""Matrix of four neighbor nodes for each node."""
shape = (ids_with_halo.shape - 2, ids_with_halo.shape - 2)
kwds = {
"strides": ids_with_halo.strides,
"buffer": ids_with_halo,
"dtype": ids_with_halo.dtype,
"offset": ids_with_halo.itemsize * (ids_with_halo.shape),
}

# kwds["offset"] = ids_with_halo.itemsize * (ids_with_halo.shape)
west_ids = np.ndarray(shape, **kwds)

kwds["offset"] = ids_with_halo.itemsize * (ids_with_halo.shape + 2)
east_ids = np.ndarray(shape, **kwds)

kwds["offset"] = ids_with_halo.itemsize
south_ids = np.ndarray(shape, **kwds)

kwds["offset"] = ids_with_halo.itemsize * (ids_with_halo.shape * 2 + 1)
north_ids = np.ndarray(shape, **kwds)

return np.vstack((east_ids.flat, north_ids.flat, west_ids.flat, south_ids.flat))

def _centered_node_ids(ids_with_halo):
"""Array of nodes taken from a matrix of nodes with a halo."""
shape = (ids_with_halo.shape - 2, ids_with_halo.shape - 2)
kwds = {
"strides": ids_with_halo.strides,
"buffer": ids_with_halo,
"dtype": ids_with_halo.dtype,
"offset": ids_with_halo.itemsize * (ids_with_halo.shape + 1),
}

return np.ndarray(shape, **kwds)

[docs]
"""Matrix of four neighbor nodes for each node."""

[docs]
):
"""Matrix of four neighbor nodes for each node."""
if open_boundary_nodes is None:
open_boundary_nodes = []

ids_with_halo = node_index_with_halo(shape, halo_indices=inactive)

# Everything that touches a closed boundary is inactive
if len(closed_boundary_nodes) > 0:
ids = _centered_node_ids(ids_with_halo)
ids.flat[closed_boundary_nodes] = inactive

neighbors = _neighbor_node_ids(ids_with_halo)

# Everything that a closed boundary touches is inactive
if len(closed_boundary_nodes) > 0:
neighbors[:, closed_boundary_nodes] = inactive

if len(open_boundary_nodes) > 0:
_set_open_boundary_neighbors(neighbors, open_boundary_nodes, inactive)

return neighbors

def _set_open_boundary_neighbors(neighbors, open_boundary_nodes, value):
"""Set values for open-boundary neighbor-nodes."""
open_boundary_neighbors = neighbors[:, open_boundary_nodes]
is_open_boundary_neighbor = _find_open_boundary_neighbors(
neighbors, open_boundary_nodes
)
nodes = np.choose(is_open_boundary_neighbor, (open_boundary_neighbors, value))
neighbors[:, open_boundary_nodes] = nodes

def _find_open_boundary_neighbors(neighbors, open_boundary_nodes):
"""Array of booleans that indicate if a neighbor is an open boundary."""
open_boundary_neighbors = neighbors[:, open_boundary_nodes]
is_open_boundary_neighbor = np.in1d(open_boundary_neighbors, open_boundary_nodes)
is_open_boundary_neighbor.shape = (neighbors.shape, len(open_boundary_nodes))
return is_open_boundary_neighbor

[docs]
def neighbor_node_array(shape, **kwds):
"""Array of neighbor nodes.

Examples
--------
>>> from landlab.utils.structured_grid import neighbor_node_array
>>> neighbors = neighbor_node_array((2, 3), inactive=-1)
>>> neighbors.T
array([[ 1,  3, -1, -1],
[ 2,  4,  0, -1],
[-1,  5,  1, -1],
[ 4, -1, -1,  0],
[ 5, -1,  3,  1],
[-1, -1,  4,  2]])
"""
closed_boundary_nodes = kwds.pop("closed_boundary_nodes", [])
open_boundary_nodes = kwds.get("open_boundary_nodes", [])

if len(closed_boundary_nodes) > 0 or len(open_boundary_nodes):
else:
neighbors = neighbor_node_ids(shape, **kwds)

return neighbors

[docs]
"""Array of neighbor cells.

Examples
--------
>>> from landlab.utils.structured_grid import neighbor_cell_array
>>> neighbors = neighbor_cell_array((2, 3), out_of_bounds=-1)
>>> len(neighbors) == 0
True

>>> neighbors = neighbor_cell_array((3, 3), out_of_bounds=-1)
>>> neighbors
array([[-1, -1, -1, -1]])

>>> neighbors = neighbor_cell_array((5, 4), out_of_bounds=-1)
>>> neighbors
array([[ 1,  2, -1, -1], [-1,  3,  0, -1],
[ 3,  4, -1,  0], [-1,  5,  2,  1],
[ 5, -1, -1,  2], [-1, -1,  4,  3]])
"""
if cell_count(shape) > 0:
shape = np.array(shape) - 2
ids = node_index_with_halo(shape, halo_indices=out_of_bounds)

neighbors = np.vstack(
(
ids[1 : shape + 1, 2:].flat,
ids[2:, 1 : shape + 1].flat,
ids[1 : shape + 1, : shape].flat,
ids[: shape, 1 : shape + 1].flat,
)
).T
if contiguous:
return neighbors.copy()
else:
return neighbors
else:
return np.array([], dtype=int)

[docs]
def diagonal_node_array(
):
"""Array of diagonal nodes.

Creates a list of IDs of the diagonal cells to each cell, as a 2D array.
Only interior cells are assigned neighbors; boundary cells get -1 for
each neighbor.  The order of the diagonal cells is [topright, topleft,
bottomleft, bottomright].

Examples
--------
>>> from landlab.utils.structured_grid import diagonal_node_array
>>> diags = diagonal_node_array((2, 3), out_of_bounds=-1)
>>> diags
array([[ 4, -1, -1, -1],
[ 5,  3, -1, -1],
[-1,  4, -1, -1],
[-1, -1, -1,  1],
[-1, -1,  0,  2],
[-1, -1,  1, -1]])
>>> diags.flags["C_CONTIGUOUS"]
True
>>> diags = diagonal_node_array((2, 3), out_of_bounds=-1, contiguous=False)
>>> diags.flags["C_CONTIGUOUS"]
False
"""
# NG didn't touch this, but she thinks this should be nodes, not cells.
ids = node_index_with_halo(shape, halo_indices=out_of_bounds)

diags = np.vstack(
(
ids[2:, 2:].flat,
ids[2:, : shape].flat,
ids[: shape, : shape].flat,
ids[: shape, 2:].flat,
)
).T

boundaries = np.empty(4, dtype=int)
diags[perimeter_nodes(shape)] = boundaries

if contiguous:
return diags.copy()
else:
return diags

[docs]
"""Array of diagonal cells.

Construct a matrix of cell indices to each diagonally adjacent cell of a
structured grid. If a cell does not have a diagonal neighbor, set the
index for that neighbor to *out_of_bounds*.

Examples
--------
A grid without any cells returns an empty array.

>>> from landlab.utils.structured_grid import diagonal_cell_array
>>> diags = diagonal_cell_array((2, 3), out_of_bounds=-1)
>>> len(diags) == 0
True

A grid that has only one cell does not have any neighbors so all of its
diagonals are set to *out_of_bounds*.
>>> diags = diagonal_cell_array((3, 3), out_of_bounds=-1)
>>> diags
array([[-1, -1, -1, -1]])

>>> diags = diagonal_cell_array((4, 4), out_of_bounds=-1)
>>> diags
array([[ 3, -1, -1, -1], [-1,  2, -1, -1],
[-1, -1, -1,  1], [-1, -1,  0, -1]])

>>> diags = diagonal_cell_array((4, 5), out_of_bounds=-1)
>>> diags
array([[ 4, -1, -1, -1], [ 5,  3, -1, -1], [-1,  4, -1, -1],
[-1, -1, -1,  1], [-1, -1,  0,  2], [-1, -1,  1, -1]])
"""
if cell_count(shape) > 0:
shape = np.array(shape) - 2
ids = node_index_with_halo(shape, halo_indices=out_of_bounds)

diags = np.vstack(
(
ids[2:, 2:].flat,
ids[2:, : shape].flat,
ids[: shape, : shape].flat,
ids[: shape, 2:].flat,
)
).T
if contiguous:
return diags.copy()
else:
return diags
else:
return np.array([], dtype=int)

[docs]
"""Array of booleans that indicate if a node has a boundary neighbor.

.. note::

DEJH thinks this method is broken since terminology update: it returns
closed neighbors, not boundary neighbors.
"""
return out_of_bounds in neighbors | out_of_bounds in diagonals

[docs]
def reshape_array(shape, array, flip_vertically=False, copy=False):
"""Reshape a flat array.

Examples
--------
>>> from landlab.utils.structured_grid import reshape_array
>>> x = np.arange(12.0)
>>> y = reshape_array((3, 4), x)
>>> y.shape == (3, 4)
True
>>> y
array([[  0.,   1.,   2.,   3.],
[  4.,   5.,   6.,   7.],
[  8.,   9.,  10.,  11.]])
>>> y.flags["C_CONTIGUOUS"]
True
>>> x = -1
>>> y[0, 0]
-1.0

>>> x = np.arange(12.0)
>>> y = reshape_array((3, 4), x, flip_vertically=True)
>>> y
array([[  8.,   9.,  10.,  11.],
[  4.,   5.,   6.,   7.],
[  0.,   1.,   2.,   3.]])
>>> y.flags["C_CONTIGUOUS"]
False
>>> x = -1
>>> y[-1, 0]
-1.0
"""
reshaped_array = array.view()

try:
reshaped_array.shape = shape
except ValueError:
raise

if flip_vertically:
flipped_array = reshaped_array[::-1, :]
if copy:
return flipped_array.copy()
else:
return flipped_array
else:
if copy:
return reshaped_array.copy()
else:
return reshaped_array

[docs]
def nodes_around_points_on_unit_grid(shape, coords, mode="raise"):
"""Array of nodes around x, y points on a grid of unit spacing.

Returns the nodes around a point on a structured grid with unit spacing
and zero origin.

Examples
--------
>>> from landlab.utils.structured_grid import nodes_around_points_on_unit_grid
>>> nodes_around_points_on_unit_grid((3, 3), (0.1, 0.1))
array([0, 3, 4, 1])

>>> nodes_around_points_on_unit_grid((3, 3), (1.0, 1.0))
array([4, 7, 8, 5])
"""
if isinstance(coords, np.ndarray):
(rows, cols) = (as_id_array(coords), as_id_array(coords))
else:
(rows, cols) = (int(coords), int(coords))

return as_id_array(
np.ravel_multi_index(
((rows, rows + 1, rows + 1, rows), (cols, cols, cols + 1, cols + 1)),
shape,
mode=mode,
).T
)

[docs]
def nodes_around_points(shape, coords, spacing=(1.0, 1.0), origin=(0.0, 0.0)):
"""Array of nodes around x, y points on a grid of non-unit spacing.

Returns the nodes around a point on a structured grid with row and column
*spacing*, and *origin*.

Examples
--------
>>> from landlab.utils.structured_grid import nodes_around_points
>>> x = np.array([0.9, 1.0])
>>> y = np.array([0.1, 1.0])
>>> nodes_around_points((3, 3), (y, x))
array([[0, 3, 4, 1],
[4, 7, 8, 5]])

>>> nodes_around_points((3, 3), (2.0, 1.0))
Traceback (most recent call last):
...
ValueError: invalid entry in coordinates array
"""
return as_id_array(
nodes_around_points_on_unit_grid(
shape,
(
(coords - origin) / spacing,
(coords - origin) / spacing,
),
)
)

[docs]
def nodes_around_point(shape, coords, spacing=(1.0, 1.0)):
"""Array of nodes around a single point on a grid of non-unit spacing."""
node_id = int(coords // spacing * shape + coords // spacing)
if node_id + shape + 1 >= shape * shape or node_id < 0:
raise ValueError("invalid entry in coordinates array")

return np.array([node_id, node_id + shape, node_id + shape + 1, node_id + 1])

```