Source code for landlab.grid.base

#! /usr/env/python
"""Python implementation of ModelGrid, a base class used to create and manage
grids for 2D numerical models.

Do NOT add new documentation here. Grid documentation is now built in a
semi-automated fashion. To modify the text seen on the web, edit the
files `docs/text_for_[gridfile].py.txt`.
"""
import fnmatch
from functools import lru_cache

import numpy as np
import xarray as xr

from landlab.utils.decorators import make_return_array_immutable

from ..core import load_params
from ..core.utils import add_module_functions_to_class
from ..field.graph_field import GraphFields
from ..layers.eventlayers import EventLayersMixIn
from ..layers.materiallayers import MaterialLayersMixIn
from ..plot.imshow import ModelGridPlotterMixIn
from ..utils.decorators import cache_result_in_object
from . import grid_funcs as gfuncs
from .decorators import (
    override_array_setitem_and_reset,
    return_id_array,
    return_readonly_id_array,
)
from .linkstatus import LinkStatus, set_status_at_link
from .nodestatus import NodeStatus

#: Indicates an index is, in some way, *bad*.
BAD_INDEX_VALUE = -1
# DEJH thinks the user should be able to override this value if they want

# Map names grid elements to the ModelGrid attribute that contains the count
# of that element in the grid.
_ARRAY_LENGTH_ATTRIBUTES = {
    "node": "number_of_nodes",
    "patch": "number_of_patches",
    "link": "number_of_links",
    "corner": "number_of_corners",
    "face": "number_of_faces",
    "cell": "number_of_cells",
    "active_link": "number_of_active_links",
    "active_face": "number_of_active_faces",
    "core_node": "number_of_core_nodes",
    "core_cell": "number_of_core_cells",
}

# Fields whose sizes can not change.
_SIZED_FIELDS = {"node", "link", "patch", "corner", "face", "cell"}


def _sort_points_into_quadrants(x, y, nodes):
    """Divide x, y points into quadrants.

    Divide points with locations given in the *x*, and *y* arrays into north,
    south, east, and west quadrants. Returns nodes contained in quadrants
    (west, east, north, south).

    Parameters
    ----------
    x : array_like
        X-coordinates of points.
    y : array_like
        Y-coordinates of points.
    nodes : array_like
        Nodes associated with points.

    Returns
    -------
    tuple of array_like
        Tuple of nodes in each coordinate. Nodes are grouped as
        (*east*, *north*, *west*, *south*).

    Examples
    --------
    >>> import numpy as np
    >>> from landlab.grid.base import _sort_points_into_quadrants
    >>> x = np.array([0, 1, 0, -1])
    >>> y = np.array([1, 0, -1, 0])
    >>> nodes = np.array([1, 2, 3, 4])
    >>> _sort_points_into_quadrants(x, y, nodes)
    (array([2]), array([1]), array([4]), array([3]))
    """
    above_x_axis = y > 0
    right_of_y_axis = x > 0
    closer_to_y_axis = np.abs(y) >= np.abs(x)

    north_nodes = nodes[above_x_axis & closer_to_y_axis]
    south_nodes = nodes[(~above_x_axis) & closer_to_y_axis]
    east_nodes = nodes[right_of_y_axis & (~closer_to_y_axis)]
    west_nodes = nodes[(~right_of_y_axis) & (~closer_to_y_axis)]

    return (east_nodes, north_nodes, west_nodes, south_nodes)


def _default_axis_names(n_dims):
    """Name of each axis.

    Parameters
    ----------
    n_dims : int
        Number of spatial dimensions.

    Returns
    -------
    tuple of str
        Name of each axis.

    Examples
    --------
    >>> from landlab.grid.base import _default_axis_names
    >>> _default_axis_names(1)
    ('x',)
    >>> _default_axis_names(2)
    ('x', 'y')
    >>> _default_axis_names(3)
    ('x', 'y', 'z')
    """
    _DEFAULT_NAMES = ("x", "y", "z")
    return _DEFAULT_NAMES[:n_dims]


def _default_axis_units(n_dims):
    """Unit names for each axis.

    Parameters
    ----------
    n_dims : int
        Number of spatial dimensions.

    Returns
    -------
    tuple of str
        Units of each axis.

    Examples
    --------
    >>> from landlab.grid.base import _default_axis_units
    >>> _default_axis_units(1)
    ('-',)
    >>> _default_axis_units(2)
    ('-', '-')
    >>> _default_axis_units(3)
    ('-', '-', '-')
    """
    return ("-",) * n_dims





[docs]class ModelGrid( GraphFields, EventLayersMixIn, MaterialLayersMixIn, ModelGridPlotterMixIn ): """Base class for 2D structured or unstructured grids for numerical models. The idea is to have at least two inherited classes, RasterModelGrid and DelaunayModelGrid, that can create and manage grids. To this might be added a GenericModelGrid, which would be an unstructured polygonal grid that doesn't necessarily obey or understand the Delaunay triangulation, but rather simply accepts an input grid from the user. Also a :class:`~.HexModelGrid` for hexagonal. Other Parameters ---------------- axis_name : tuple, optional Name of axes axis_units : tuple, optional Units of coordinates """ #: Indicates a node is *bad index*. BAD_INDEX = BAD_INDEX_VALUE #: Indicates a node is *core*. BC_NODE_IS_CORE = NodeStatus.CORE #: Indicates a boundary node has a fixed value. BC_NODE_IS_FIXED_VALUE = NodeStatus.FIXED_VALUE #: Indicates a boundary node has a fixed gradient. BC_NODE_IS_FIXED_GRADIENT = NodeStatus.FIXED_GRADIENT #: Indicates a boundary node is wrap-around. BC_NODE_IS_LOOPED = NodeStatus.LOOPED #: Indicates a boundary node is closed BC_NODE_IS_CLOSED = NodeStatus.CLOSED #: Indicates a link is *active*, and can carry flux BC_LINK_IS_ACTIVE = LinkStatus.ACTIVE #: Indicates a link has a fixed gradient value, and behaves as a boundary BC_LINK_IS_FIXED = LinkStatus.FIXED #: Indicates a link is *inactive*, and cannot carry flux BC_LINK_IS_INACTIVE = LinkStatus.INACTIVE #: Grid elements on which fields can be placed. VALID_LOCATIONS = ("node", "link", "patch", "corner", "face", "cell", "grid") at_node = {} # : Values defined at nodes at_link = {} # : Values defined at links at_patch = {} # : Values defined at patches at_corner = {} # : Values defined at corners at_face = {} # : Values defined at faces at_cell = {} # : Values defined at cells at_grid = {} # : Values defined at grid
[docs] @classmethod def from_file(cls, file_like): """Create grid from a file-like object. File to load either as a file-like object, path to an existing file, or the contents of a file as a string. Parameters ---------- file_like : File-like object, filepath, or string. Examples -------- >>> from io import StringIO >>> from landlab import RasterModelGrid >>> filelike = StringIO(''' ... shape: ... - 3 ... - 4 ... xy_spacing: 2 ... ''') >>> grid = RasterModelGrid.from_file(filelike) >>> grid.x_of_node array([ 0., 2., 4., 6., 0., 2., 4., 6., 0., 2., 4., 6.]) >>> grid.y_of_node array([ 0., 0., 0., 0., 2., 2., 2., 2., 4., 4., 4., 4.]) """ params = load_params(file_like) return cls.from_dict(params)
[docs] @classmethod def from_dict(cls, params): """Create grid from dictionary. Parameters ---------- params : dictionary Dictionary of required parameters to create a model grid. Examples -------- >>> from landlab import RasterModelGrid >>> params = {"shape": (3,4), "xy_spacing": 2} >>> grid = RasterModelGrid.from_dict(params) >>> grid.x_of_node array([ 0., 2., 4., 6., 0., 2., 4., 6., 0., 2., 4., 6.]) >>> grid.y_of_node array([ 0., 0., 0., 0., 2., 2., 2., 2., 4., 4., 4., 4.]) """ return cls(**params)
def __init__(self, **kwds): axis_units = kwds.pop("xy_axis_units", "-") axis_name = kwds.pop("xy_axis_name", ("x", "y")) super().__init__() self.new_field_location("node", self.number_of_nodes) self.new_field_location("link", self.number_of_links) self.new_field_location("patch", self.number_of_patches) self.new_field_location("corner", self.number_of_corners) self.new_field_location("face", self.number_of_faces) self.new_field_location("cell", self.number_of_cells) self.new_field_location("grid", None) self.default_group = "node" # self.axis_name = kwds.get("axis_name", _default_axis_names(self.ndim)) # self.axis_units = kwds.get("axis_units", _default_axis_units(self.ndim)) self._ref_coord = tuple(kwds.get("xy_of_reference", (0.0, 0.0))) self._link_length = None self._all_node_distances_map = None self._all_node_azimuths_map = None self.bc_set_code = 0 self._axis_units = tuple(np.broadcast_to(axis_units, self.ndim)) self._axis_name = tuple(np.broadcast_to(axis_name, self.ndim))
[docs] def fields(self, include="*", exclude=None): """List of fields held by the grid. The returned field names are returned as their canonical names. This is, as a string of the for "at_<location>:<name>". This allows for fields with the same name to be defined at different grid locations. You could have, for example, a variable "elevation" defined at both *nodes* and *links*. Both the *include* and *exclude* patterns are glob-style expressions, not regular expressions. If either *include* or *exclude* are lists, then the patterns are matched using an "or". The *include* filters are applied before the *exclude* filters. Parameters ---------- include : str, or iterable of str, optional Glob-style pattern for field names to include. exclude : str, or iterable of str, optional Glob-style pattern for field names to exclude. Returns ------- set Filtered set of canonical field names held by the grid Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((3, 4)) Add some fields to the grid. Notice that we have defined a field named "elevation" at both *nodes* and *links*. >>> _ = grid.add_full("elevation", 3.0, at="node") >>> _ = grid.add_full("elevation", 4.0, at="link") >>> _ = grid.add_full("temperature", 5.0, at="node") >>> sorted(grid.fields()) ['at_link:elevation', 'at_node:elevation', 'at_node:temperature'] >>> sorted(grid.fields(include="at_node*")) ['at_node:elevation', 'at_node:temperature'] >>> sorted(grid.fields(include="at_node*", exclude="*temp*")) ['at_node:elevation'] Fields can also be defined at *layers*. In the following example we've filtered the results to just return the layer fields. >>> grid.event_layers.add(1.0, rho=0.5) >>> sorted(grid.fields(include="at_layer*")) ['at_layer:rho'] If a list, the fields are matched with an "or". >>> sorted(grid.fields(include=["at_node*", "*elevation*"])) ['at_link:elevation', 'at_node:elevation', 'at_node:temperature'] """ if isinstance(include, str): include = [include] if isinstance(exclude, str): exclude = [exclude] layer_groups = set(["_".join(["layer", at]) for at in self.groups]) layer_groups.add("layer") canonical_names = set() for at in self.groups | layer_groups: try: canonical_names.update( ["at_{0}:{1}".format(at, name) for name in self[at]] ) except KeyError: pass names = set() for pattern in include: names.update(fnmatch.filter(canonical_names, pattern)) for pattern in exclude or []: names.difference_update(fnmatch.filter(canonical_names, pattern)) return names
[docs] def as_dataarray(self, name, at=None, time=None): """Create an xarray DataArray representation of a grid field. Parameters ---------- name : str Name of a field. This can either be a canonical field name (of the form "at_<element>:<field_name>", or just the field name. In the latter case, use the *at* keyword to specify where the field is defined. at : str, optional The grid elements on which the field is defined. Use this only if *name* is not a canonical field name that already contains the grid element. Returns ------- DataArray The field represented as a newly-created *xarray* *DataArray*. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((3, 4)) >>> _ = grid.add_full("elevation", 3.0, at="node") >>> grid.as_dataarray("at_node:elevation") <xarray.DataArray 'at_node:elevation' (node: 12)> array([ 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3., 3.]) Dimensions without coordinates: node >>> all( ... grid.as_dataarray("at_node:elevation") ... == grid.as_dataarray("elevation", at="node") ... ) True """ if at is None: at, field_name = name[len("at_") :].split(":") else: field_name = name values = getattr(self, f"at_{at}")[field_name] if at == "layer": # For backward compatibility at = "layer_cell" dims = at.split("_") if at == "grid": dims = [f"{field_name}_per_grid"] if values.shape else [] if time is not None and "layer" not in dims: dims = ["time"] + dims values = values[None] return xr.DataArray(values, dims=dims, name=f"at_{at}:{field_name}")
[docs] def as_dataset(self, include="*", exclude=None, time=None): """Create an xarray Dataset representation of a grid. This method creates a new xarray Dataset object that contains the grid's data fields. A particular grid type (e.g. a *RasterModelGrid*) should define its own *as_dataset* method to represent that particular grid type as a *Dataset* and then call this method to add the data fields. Parameters ---------- include : str or iterable or str Glob-style patterns of fields to include in the dataset. exclude : str or iterable or str Glob-style patterns of fields to exclude from the dataset. Returns ------- Dataset An xarray Dataset representation of a *ModelGrid*. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((3, 4)) Add some fields to the grid. Notice that we have defined a field named "elevation" at both *nodes* and *links*. >>> _ = grid.add_full("elevation", 3.0, at="node") >>> _ = grid.add_full("elevation", 4.0, at="link") >>> _ = grid.add_full("temperature", 5.0, at="node") >>> ds = grid.as_dataset() >>> sorted(ds.dims.items()) [('dim', 2), ('link', 17), ('node', 12)] >>> sorted([var for var in ds.data_vars if var.startswith("at_")]) ['at_link:elevation', 'at_node:elevation', 'at_node:temperature'] >>> grid.event_layers.add(1.0, rho=0.5) >>> ds = grid.as_dataset() >>> sorted(ds.dims.items()) [('cell', 2), ('dim', 2), ('layer', 1), ('link', 17), ('node', 12)] >>> sorted([var for var in ds.data_vars if var.startswith("at_")]) ['at_layer_cell:rho', 'at_layer_cell:thickness', 'at_link:elevation', 'at_node:elevation', 'at_node:temperature'] """ names = self.fields(include=include, exclude=exclude) data = {} for name in names: array = self.as_dataarray(name, time=time) data[array.name] = array if self.event_layers.number_of_layers > 0: data["at_layer_cell:thickness"] = (("layer", "cell"), self.event_layers.dz) data["status_at_node"] = (("node",), self.status_at_node) if time is not None: data["time"] = (("time",), [time]) return xr.Dataset(data)
@property def xy_of_reference(self): """Return the coordinates (x, y) of the reference point. For RasterModelGrid and HexModelGrid the reference point is the minimum of x_of_node and of y_of_node. By default it is (0, 0). For VoronoiDelaunayGrid the reference point is (0, 0). For RadialModelGrid it is the (x, y) of the center point. The intention of these coordinates is to provide a method to store the large float values of projected coordinates. Example ------- >>> from landlab import RasterModelGrid >>> rmg = RasterModelGrid((4, 5), ... xy_of_reference = (12345, 678910)) >>> rmg.xy_of_reference (12345, 678910) >>> rmg.xy_of_reference = (98765, 43210) >>> rmg.xy_of_reference (98765, 43210) """ return self._ref_coord @xy_of_reference.setter def xy_of_reference(self, new_xy_of_reference): """Set a new value for the model grid xy_of_reference.""" self._ref_coord = (new_xy_of_reference[0], new_xy_of_reference[1]) @property def ndim(self): """Number of spatial dimensions of the grid. LLCATS: GINF """ return 2 @property @override_array_setitem_and_reset("reset_status_at_node") def status_at_node(self): """Get array of the boundary status for each node. Examples -------- >>> import numpy as np >>> from landlab import LinkStatus, NodeStatus, RasterModelGrid >>> mg = RasterModelGrid((4, 5)) >>> mg.status_at_node.reshape((4, 5)) array([[1, 1, 1, 1, 1], [1, 0, 0, 0, 1], [1, 0, 0, 0, 1], [1, 1, 1, 1, 1]], dtype=uint8) >>> np.any(mg.status_at_link == LinkStatus.FIXED) False >>> mg.status_at_node[mg.nodes_at_left_edge] = NodeStatus.FIXED_GRADIENT >>> mg.status_at_node.reshape((4, 5)) array([[2, 1, 1, 1, 1], [2, 0, 0, 0, 1], [2, 0, 0, 0, 1], [2, 1, 1, 1, 1]], dtype=uint8) >>> np.any(mg.status_at_link == LinkStatus.FIXED) # links auto-update True LLCATS: NINF BC """ return self._node_status @status_at_node.setter def status_at_node(self, new_status): """Set the array of node boundary statuses.""" self._node_status[:] = new_status[:] self.reset_status_at_node() @property @cache_result_in_object() @return_readonly_id_array def active_adjacent_nodes_at_node(self): """Adjacent nodes for each grid node. For each grid node, get the adjacent nodes ordered counterclockwise starting from the positive x axis. Examples -------- >>> from landlab import RasterModelGrid, HexModelGrid >>> grid = RasterModelGrid((4, 5)) >>> grid.active_adjacent_nodes_at_node[(-1, 6, 2), ] array([[-1, -1, -1, -1], [ 7, 11, 5, 1], [-1, 7, -1, -1]]) Setting a node to closed causes all links touching it to be inactive. >>> grid.status_at_node[6] = grid.BC_NODE_IS_CLOSED >>> grid.active_adjacent_nodes_at_node[(-1, 6, 2), ] array([[-1, -1, -1, -1], [-1, -1, -1, -1], [-1, 7, -1, -1]]) >>> grid.active_adjacent_nodes_at_node[7] array([ 8, 12, -1, 2]) >>> grid.active_adjacent_nodes_at_node[2] array([-1, 7, -1, -1]) >>> grid = HexModelGrid((3, 2)) >>> grid.status_at_node[0] = grid.BC_NODE_IS_CLOSED >>> grid.active_adjacent_nodes_at_node array([[-1, -1, -1, -1, -1, -1], [-1, 3, -1, -1, -1, -1], [ 3, -1, -1, -1, -1, -1], [ 4, 6, 5, 2, -1, 1], [-1, 3, -1, -1, -1, -1], [-1, -1, 3, -1, -1, -1], [-1, 3, -1, -1, -1, -1]]) LLCATS: NINF CONN BC """ return np.choose( self.status_at_link[self.links_at_node] == LinkStatus.ACTIVE, (-1, self.adjacent_nodes_at_node), ) @property @make_return_array_immutable @cache_result_in_object() def active_link_dirs_at_node(self): """Link flux directions at each node: 1=incoming flux, -1=outgoing flux, 0=no flux. Note that inactive links receive zero, but active and fixed links are both reported normally. Returns ------- (NODES, LINKS) ndarray of int Link directions relative to the nodes of a grid. The shape of the matrix will be number of nodes rows by max number of links per node. A zero indicates no link at this position. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((4, 3)) >>> grid.status_at_node[grid.nodes_at_left_edge] = grid.BC_NODE_IS_CLOSED >>> grid.active_link_dirs_at_node # doctest: +NORMALIZE_WHITESPACE array([[ 0, 0, 0, 0], [ 0, -1, 0, 0], [ 0, 0, 0, 0], [ 0, 0, 0, 0], [-1, -1, 0, 1], [ 0, 0, 1, 0], [ 0, 0, 0, 0], [-1, -1, 0, 1], [ 0, 0, 1, 0], [ 0, 0, 0, 0], [ 0, 0, 0, 1], [ 0, 0, 0, 0]], dtype=int8) LLCATS: NINF LINF CONN """ return np.choose( self.link_status_at_node == LinkStatus.ACTIVE, (0, self.link_dirs_at_node) ) @property @make_return_array_immutable @cache_result_in_object() def link_status_at_node(self): return self.status_at_link[self.links_at_node] @property @return_readonly_id_array @cache_result_in_object() def core_nodes(self): """Get array of core nodes. Examples -------- >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((4, 5)) >>> mg.core_nodes array([ 6, 7, 8, 11, 12, 13]) LLCATS: NINF BC """ return np.where(self.status_at_node == NodeStatus.CORE)[0] @property @return_readonly_id_array def boundary_nodes(self): """Get array of boundary nodes. Examples -------- >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((4, 5)) >>> mg.boundary_nodes array([ 0, 1, 2, 3, 4, 5, 9, 10, 14, 15, 16, 17, 18, 19]) LLCATS: NINF BC """ try: return self._boundary_nodes except AttributeError: (boundary_node_ids,) = np.where(self._node_status != NodeStatus.CORE) return boundary_node_ids @property @return_readonly_id_array def open_boundary_nodes(self): """Get array of open boundary nodes. Examples -------- >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((4, 5)) >>> for edge in (mg.nodes_at_left_edge, mg.nodes_at_right_edge, ... mg.nodes_at_bottom_edge): ... mg.status_at_node[edge] = mg.BC_NODE_IS_CLOSED >>> mg.open_boundary_nodes array([16, 17, 18]) LLCATS: NINF BC """ (open_boundary_node_ids,) = np.where( (self._node_status != self.BC_NODE_IS_CLOSED) & (self._node_status != self.BC_NODE_IS_CORE) ) return open_boundary_node_ids @property @return_readonly_id_array def closed_boundary_nodes(self): """Get array of closed boundary nodes. Examples -------- >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((4, 5)) >>> mg.status_at_node[mg.nodes_at_top_edge] = mg.BC_NODE_IS_CLOSED >>> mg.closed_boundary_nodes array([15, 16, 17, 18, 19]) LLCATS: NINF BC """ (closed_boundary_node_ids,) = np.where( self._node_status == self.BC_NODE_IS_CLOSED ) return closed_boundary_node_ids @property @return_readonly_id_array def fixed_gradient_boundary_nodes(self): """Get array of fixed gradient boundary nodes. Examples -------- >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((4, 5)) >>> mg.status_at_node[mg.nodes_at_top_edge] = mg.BC_NODE_IS_FIXED_GRADIENT >>> mg.fixed_gradient_boundary_nodes array([15, 16, 17, 18, 19]) LLCATS: NINF BC """ (fixed_gradient_boundary_node_ids,) = np.where( self._node_status == self.BC_NODE_IS_FIXED_GRADIENT ) return fixed_gradient_boundary_node_ids @property @return_readonly_id_array def fixed_gradient_boundary_node_fixed_link(self): """An array of the fixed_links connected to fixed gradient boundary nodes. Note that on a raster, some nodes (notably the corners) can be `NodeStatus.FIXED_GRADIENT`, but not have a true `LinkStatus.FIXED` neighboring link. In such cases, the link returned will be a closed link joining the corner node to a neighboring `NodeStatus.FIXED_GRADIENT` node (see example). An AssertionError will be raised if for some reason a `NodeStatus.FIXED_GRADIENT` node exists which has neither a `NodeStatus.FIXED_GRADIENT` neighbor, or a `LinkStatus.FIXED`. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((3, 4)) >>> leftedge = grid.nodes_at_left_edge >>> grid.status_at_node[leftedge] = grid.BC_NODE_IS_FIXED_GRADIENT >>> grid.fixed_gradient_boundary_nodes array([0, 4, 8]) >>> grid.fixed_gradient_boundary_node_fixed_link array([ 3, 7, 10]) """ try: return self._fixed_gradient_boundary_node_links except AttributeError: self._create_fixed_gradient_boundary_node_links() return self._fixed_gradient_boundary_node_links @property @return_readonly_id_array def fixed_gradient_boundary_node_anchor_node(self): """Returns the node at the other end of the fixed link for a fixed gradient boundary node. Degenerate `NodeStatus.FIXED_GRADIENT` nodes (e.g., corners) are handled as in :func:`fixed_gradient_boundary_node_fixed_link`, by pointing to a neighboring `NodeStatus.FIXED_GRADIENT` node. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((3, 4)) >>> leftedge = grid.nodes_at_left_edge >>> grid.status_at_node[leftedge] = grid.BC_NODE_IS_FIXED_GRADIENT >>> grid.fixed_gradient_boundary_nodes array([0, 4, 8]) >>> grid.fixed_gradient_boundary_node_fixed_link array([ 3, 7, 10]) >>> grid.fixed_gradient_boundary_node_anchor_node array([4, 5, 4]) """ try: return self._fixed_gradient_boundary_node_anchor_node except AttributeError: self._create_fixed_gradient_boundary_node_anchor_node() return self._fixed_gradient_boundary_node_anchor_node def _create_fixed_gradient_boundary_node_links(self): """Builds a data structure to hold the fixed_links which control the values of any `NodeStatus.FIXED_GRADIENT` nodes in the grid. An AssertionError will be raised if for some reason a `NodeStatus.FIXED_GRADIENT` node exists which has neither a `NodeStatus.FIXED_GRADIENT` neighbor, or a `LinksStatus.FIXED`. """ self._fixed_grad_links_created = True self._fixed_gradient_boundary_node_links = np.empty_like( self.fixed_gradient_boundary_nodes, dtype=int ) fix_nodes = self.fixed_gradient_boundary_nodes neighbor_links = self.links_at_node[fix_nodes] # -1s boundary_exists = self.link_dirs_at_node[fix_nodes] # next line retains -1 indexes link_stat_badind = self.status_at_link[neighbor_links] == LinkStatus.FIXED true_connection = np.logical_and(link_stat_badind, boundary_exists) true_fix_nodes = true_connection.sum(axis=1).astype(bool) self._fixed_gradient_boundary_node_links[true_fix_nodes] = neighbor_links[ true_connection ] # resolve any corner nodes neighbor_nodes = self.adjacent_nodes_at_node[fix_nodes] # BAD_INDEX_VALUEs neighbor_nodes[neighbor_nodes == self.BAD_INDEX] = -1 fixed_grad_neighbor = np.logical_and( (self.status_at_node[neighbor_nodes] == NodeStatus.FIXED_GRADIENT), boundary_exists, ) # ^True when NodeStatus.FIXED_GRADIENT for real # winnow it down to only one possibility for fixed_grad neighbor: which_neighbor = np.argmax(fixed_grad_neighbor, axis=1) indexing_range = np.arange(fixed_grad_neighbor.shape[0]) a_link_to_fixed_grad = neighbor_links[indexing_range, which_neighbor] corners = np.logical_not(true_fix_nodes) assert np.all(fixed_grad_neighbor[indexing_range, which_neighbor][corners]) self._fixed_gradient_boundary_node_links[corners] = a_link_to_fixed_grad[ corners ] def _create_fixed_gradient_boundary_node_anchor_node(self): """Builds a data structure to hold the nodes which anchor the values of any `NodeStatus.FIXED_GRADIENT` nodes in the grid, i.e., those at the other ends of the `LinkStatus.FIXED`. An AssertionError will be raised if for some reason a `NodeStatus.FIXED_GRADIENT` node exists which has neither a `NodeStatus.FIXED_GRADIENT` neighbor, or a `LinkStatus.FIXED`. """ self._fixed_grad_links_created = True fix_grad_nodes = self.fixed_gradient_boundary_nodes self._fixed_gradient_boundary_node_anchor_node = np.empty_like(fix_grad_nodes) heads_and_tails = np.empty((fix_grad_nodes.size, 2)) which_one = np.empty_like(heads_and_tails, dtype=bool) heads_and_tails[:, 0] = self.node_at_link_head[ self.fixed_gradient_boundary_node_fixed_link ] heads_and_tails[:, 1] = self.node_at_link_tail[ self.fixed_gradient_boundary_node_fixed_link ] which_one[:, 0] = heads_and_tails[:, 0] == fix_grad_nodes which_one[:, 1] = heads_and_tails[:, 1] == fix_grad_nodes assert np.all(which_one.sum(axis=1) == 1) self._fixed_gradient_boundary_node_anchor_node = heads_and_tails[ np.logical_not(which_one) ] @property @return_readonly_id_array @cache_result_in_object() def fixed_value_boundary_nodes(self): """Get array of fixed value boundary nodes. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((4, 5)) Initially all the perimeter nodes are fixed value boundary. >>> grid.fixed_value_boundary_nodes array([ 0, 1, 2, 3, 4, 5, 9, 10, 14, 15, 16, 17, 18, 19]) Set left, right, and bottom edges to closed. >>> for edge in (grid.nodes_at_left_edge, grid.nodes_at_right_edge, ... grid.nodes_at_bottom_edge): ... grid.status_at_node[edge] = grid.BC_NODE_IS_CLOSED Now nodes on just the top edge are fixed. >>> grid.fixed_value_boundary_nodes array([16, 17, 18]) LLCATS: NINF BC """ return np.where(self._node_status == NodeStatus.FIXED_VALUE)[0] @property @return_readonly_id_array @cache_result_in_object() def active_faces(self): """Get array of active faces. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((3, 4)) >>> grid.active_faces array([0, 1, 2, 3, 4, 5, 6]) >>> grid.status_at_node[6] = grid.BC_NODE_IS_CLOSED >>> grid.active_faces array([0, 2, 5]) LLCATS: FINF BC """ return self.face_at_link[self.active_links] @property @return_readonly_id_array @cache_result_in_object() def active_links(self): """Get array of active links. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((3, 4)) >>> grid.active_links array([ 4, 5, 7, 8, 9, 11, 12]) LLCATS: LINF BC """ return np.where(self.status_at_link == LinkStatus.ACTIVE)[0] @property @return_readonly_id_array @cache_result_in_object() def fixed_links(self): """Get array of fixed links. Examples -------- >>> from landlab import NodeStatus, RasterModelGrid >>> grid = RasterModelGrid((3, 4)) >>> grid.status_at_node # doctest: +NORMALIZE_WHITESPACE array([1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1], dtype=uint8) >>> grid.fixed_links.size 0 >>> grid.status_at_node[:4] = NodeStatus.FIXED_GRADIENT >>> grid.status_at_node # doctest: +NORMALIZE_WHITESPACE array([2, 2, 2, 2, 1, 0, 0, 1, 1, 1, 1, 1], dtype=uint8) >>> grid.fixed_links array([4, 5]) LLCATS: LINF BC """ return np.where(self.status_at_link == LinkStatus.FIXED)[0] @property @cache_result_in_object() @return_readonly_id_array def node_at_core_cell(self): """Get array of nodes associated with core cells. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((4, 5)) Initially each cell's node is core. >>> grid.node_at_core_cell array([ 6, 7, 8, 11, 12, 13]) Setting a node to closed causes means its cell is also "closed". >>> grid.status_at_node[8] = grid.BC_NODE_IS_CLOSED >>> grid.node_at_core_cell array([ 6, 7, 11, 12, 13]) LLCATS: NINF CINF BC CONN """ return np.where(self.status_at_node == NodeStatus.CORE)[0] @property @make_return_array_immutable @cache_result_in_object() def core_cells(self): """Get array of core cells. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((4, 5)) Initially all of the cells are "core". >>> grid.core_cells array([0, 1, 2, 3, 4, 5]) Setting a node to closed causes its cell to no longer be core. >>> grid.status_at_node[8] = grid.BC_NODE_IS_CLOSED >>> grid.core_cells array([0, 1, 3, 4, 5]) LLCATS: CINF BC """ return self.cell_at_node[self.core_nodes] @property def number_of_active_faces(self): """Total number of active faces. Returns ------- int Total number of active faces in the grid. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((3, 4)) >>> grid.number_of_active_faces 7 The number of active faces is updated when a node status changes. >>> grid.status_at_node[6] = grid.BC_NODE_IS_CLOSED >>> grid.number_of_active_faces 3 LLCATS: FINF BC """ return self.active_faces.size @property def number_of_core_nodes(self): """Number of core nodes. The number of core nodes on the grid (i.e., excluding all boundary nodes). Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((4, 5)) >>> grid.number_of_core_nodes 6 >>> grid.status_at_node[7] = grid.BC_NODE_IS_CLOSED >>> grid.number_of_core_nodes 5 LLCATS: NINF BC """ return self.core_nodes.size @property def number_of_core_cells(self): """Number of core cells. A core cell excludes all boundary cells. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((4, 5)) >>> grid.number_of_core_cells 6 >>> grid.status_at_node[7] = grid.BC_NODE_IS_CLOSED >>> grid.number_of_core_cells 5 LLCATS: CINF BC """ return self.core_cells.size @property def number_of_active_links(self): """Number of active links. Examples -------- >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((4, 5)) >>> mg.number_of_active_links 17 >>> for edge in (mg.nodes_at_left_edge, mg.nodes_at_right_edge, ... mg.nodes_at_bottom_edge): ... mg.status_at_node[edge] = mg.BC_NODE_IS_CLOSED >>> mg.number_of_active_links 10 LLCATS: LINF BC """ return self.active_links.size @property def number_of_fixed_links(self): """Number of fixed links. Examples -------- >>> from landlab import NodeStatus, RasterModelGrid >>> mg = RasterModelGrid((4, 5)) >>> mg.number_of_fixed_links 0 >>> mg.status_at_node[mg.nodes_at_top_edge] = NodeStatus.FIXED_GRADIENT >>> mg.number_of_fixed_links 3 LLCATS: LINF BC """ return self.fixed_links.size
[docs] def number_of_elements(self, name): """Number of instances of an element. Get the number of instances of a grid element in a grid. Parameters ---------- name : {'node', 'cell', 'link', 'face', 'core_node', 'core_cell', 'active_link', 'active_face'} Name of the grid element. Returns ------- int Number of elements in the grid. Examples -------- >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((4, 5)) >>> mg.number_of_elements('node') 20 >>> mg.number_of_elements('core_cell') 6 >>> mg.number_of_elements('link') 31 >>> mg.number_of_elements('active_link') 17 >>> mg.status_at_node[8] = mg.BC_NODE_IS_CLOSED >>> mg.number_of_elements('link') 31 >>> mg.number_of_elements('active_link') 13 LLCATS: GINF """ try: return getattr(self, _ARRAY_LENGTH_ATTRIBUTES[name]) except KeyError: raise TypeError("{name}: element name not understood".format(name=name))
[docs] @make_return_array_immutable def node_axis_coordinates(self, axis=0): """Get the coordinates of nodes along a particular axis. Return node coordinates from a given *axis* (defaulting to 0). Axis numbering is the same as that for numpy arrays. That is, the zeroth axis is along the rows, and the first along the columns. Parameters ---------- axis : int, optional Coordinate axis. Returns ------- ndarray Coordinates of nodes for a given axis. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((4, 5)) >>> grid.node_axis_coordinates(0) # doctest: +NORMALIZE_WHITESPACE array([ 0., 0., 0., 0., 0., 1., 1., 1., 1., 1., 2., 2., 2., 2., 2., 3., 3., 3., 3., 3.]) >>> grid.node_axis_coordinates(1) # doctest: +NORMALIZE_WHITESPACE array([ 0., 1., 2., 3., 4., 0., 1., 2., 3., 4., 0., 1., 2., 3., 4., 0., 1., 2., 3., 4.]) LLCATS: GINF NINF MEAS """ AXES = ("node_y", "node_x") try: return getattr(self, AXES[axis]) except IndexError: raise ValueError("'axis' entry is out of bounds")
@property def axis_units(self): """Get units for each axis. Returns ------- tuple of str The units (as a string) for each of a grid's coordinates. Examples -------- >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((4, 5), xy_spacing=(3., 2.)) >>> mg.axis_units ('-', '-') >>> mg.axis_units = ("degrees_north", "degrees_east") >>> mg.axis_units ('degrees_north', 'degrees_east') >>> mg.axis_units = "m" >>> mg.axis_units ('m', 'm') LLCATS: GINF """ return self._axis_units @axis_units.setter def axis_units(self, new_units): """Set the units for each coordinate axis.""" new_units = np.broadcast_to(new_units, self.ndim) if len(new_units) != self.ndim: raise ValueError("length of units does not match grid dimension") self._axis_units = tuple(new_units) @property def axis_name(self): """Get the name of each coordinate axis. Returns ------- tuple of str The names of each axis. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((4, 5)) >>> grid.axis_name ('x', 'y') >>> grid.axis_name = ('lon', 'lat') >>> grid.axis_name ('lon', 'lat') LLCATS: GINF """ return self._axis_name @axis_name.setter def axis_name(self, new_names): """Set the names of a grid's coordinate axes. Raises ------ ValueError If the number of dimension do not match. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((4, 5)) >>> grid.axis_name = ("y", "x") >>> grid.axis_name = ("lon", "lat") >>> grid.axis_name ('lon', 'lat') """ new_names = np.broadcast_to(new_names, self.ndim) if len(new_names) != self.ndim: raise ValueError("length of names does not match grid dimension") self._axis_name = tuple(new_names) @property @make_return_array_immutable @cache_result_in_object() def status_at_link(self): """Get array of the status of all links. Examples -------- >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((4, 5)) >>> mg.status_at_node[mg.nodes_at_left_edge] = mg.BC_NODE_IS_CLOSED >>> mg.status_at_node[mg.nodes_at_right_edge] = mg.BC_NODE_IS_FIXED_GRADIENT >>> mg.status_at_link # doctest: +NORMALIZE_WHITESPACE array([4, 4, 4, 4, 4, 0, 0, 0, 4, 4, 0, 0, 2, 4, 0, 0, 0, 4, 4, 0, 0, 2, 4, 0, 0, 0, 4, 4, 4, 4, 4], dtype=uint8) LLCATS: BC LINF """ return set_status_at_link(self.status_at_node[self.nodes_at_link]) @property @make_return_array_immutable @cache_result_in_object() def angle_of_link_about_head(self): """Find and return the angle of a link about the node at the link head. Because links have direction, their angle can be specified as an angle about either the node at the link head, or the node at the link tail. The default behaviour of `angle_of_link` is to return the angle about the link tail, but this method gives the angle about the link head. Examples -------- >>> from landlab import HexModelGrid >>> import numpy as np >>> grid = HexModelGrid((3, 2), node_layout="hex") >>> np.round(grid.angle_of_link[:3] / np.pi * 3.0) array([ 0., 2., 1.]) >>> np.round(grid.angle_of_link_about_head[:3] / np.pi * 3.0) # 60 deg segments array([ 3., 5., 4.]) LLCATS: LINF MEAS """ angles = np.arctan2(-np.sin(self.angle_of_link), -np.cos(self.angle_of_link)) return np.mod(angles, 2.0 * np.pi, out=angles) @property @make_return_array_immutable def patches_present_at_node(self): """A boolean array, False where a patch has a closed node or is missing. The array is the same shape as :func:`patches_at_node`, and is designed to mask it. Note that in cases where patches may have more than 3 nodes (e.g., rasters), a patch is considered still present as long as at least 3 open nodes are present. Examples -------- >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((3, 3)) >>> mg.status_at_node[mg.nodes_at_top_edge] = mg.BC_NODE_IS_CLOSED >>> mg.patches_at_node array([[ 0, -1, -1, -1], [ 1, 0, -1, -1], [-1, 1, -1, -1], [ 2, -1, -1, 0], [ 3, 2, 0, 1], [-1, 3, 1, -1], [-1, -1, -1, 2], [-1, -1, 2, 3], [-1, -1, 3, -1]]) >>> mg.patches_present_at_node array([[ True, False, False, False], [ True, True, False, False], [False, True, False, False], [False, False, False, True], [False, False, True, True], [False, False, True, False], [False, False, False, False], [False, False, False, False], [False, False, False, False]], dtype=bool) >>> 1 in mg.patches_at_node * mg.patches_present_at_node True >>> 2 in mg.patches_at_node * mg.patches_present_at_node False LLCATS: PINF NINF """ try: return self._patches_present_mask except AttributeError: self.patches_at_node self._reset_patch_status() return self._patches_present_mask @property @make_return_array_immutable def patches_present_at_link(self): """A boolean array, False where a patch has a closed node or is missing. The array is the same shape as :func:`patches_at_link`, and is designed to mask it. Examples -------- >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((3, 3)) >>> mg.status_at_node[mg.nodes_at_top_edge] = mg.BC_NODE_IS_CLOSED >>> mg.patches_at_link array([[-1, 0], [-1, 1], [ 0, -1], [ 1, 0], [-1, 1], [ 0, 2], [ 1, 3], [ 2, -1], [ 3, 2], [-1, 3], [ 2, -1], [ 3, -1]]) >>> mg.patches_present_at_link array([[False, True], [False, True], [ True, False], [ True, True], [False, True], [ True, False], [ True, False], [False, False], [False, False], [False, False], [False, False], [False, False]], dtype=bool) >>> 1 in mg.patches_at_link * mg.patches_present_at_link True >>> 2 in mg.patches_at_link * mg.patches_present_at_link False LLCATS: PINF LINF """ try: return self._patches_present_link_mask except AttributeError: self.patches_at_node self._reset_patch_status() return self._patches_present_link_mask @property @make_return_array_immutable def number_of_patches_present_at_node(self): """Return the number of patches at a node without a closed node. Examples -------- >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((3, 3)) >>> mg.status_at_node[mg.nodes_at_top_edge] = mg.BC_NODE_IS_CLOSED >>> mg.patches_present_at_node array([[ True, False, False, False], [ True, True, False, False], [False, True, False, False], [False, False, False, True], [False, False, True, True], [False, False, True, False], [False, False, False, False], [False, False, False, False], [False, False, False, False]], dtype=bool) >>> mg.number_of_patches_present_at_node array([1, 2, 1, 1, 2, 1, 0, 0, 0]) LLCATS: PINF NINF BC """ try: return self._number_of_patches_present_at_node except AttributeError: self.patches_at_node self._reset_patch_status() return self._number_of_patches_present_at_node @property @make_return_array_immutable def number_of_patches_present_at_link(self): """Return the number of patches at a link without a closed node. Examples -------- >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((3, 3)) >>> mg.status_at_node[mg.nodes_at_top_edge] = mg.BC_NODE_IS_CLOSED >>> mg.patches_present_at_link array([[False, True], [False, True], [ True, False], [ True, True], [False, True], [ True, False], [ True, False], [False, False], [False, False], [False, False], [False, False], [False, False]], dtype=bool) >>> mg.number_of_patches_present_at_link array([1, 1, 1, 2, 1, 1, 1, 0, 0, 0, 0, 0]) LLCATS: PINF LINF BC """ try: return self._number_of_patches_present_at_link except AttributeError: self.patches_at_node self._reset_patch_status() return self._number_of_patches_present_at_link def _reset_patch_status(self): """Creates the array which stores patches_present_at_node. Call whenever boundary conditions are updated on the grid. """ from landlab import RasterModelGrid, VoronoiDelaunayGrid node_status_at_patch = self.status_at_node[self.nodes_at_patch] if isinstance(self, RasterModelGrid): max_nodes_at_patch = 4 elif isinstance(self, VoronoiDelaunayGrid): max_nodes_at_patch = 3 else: max_nodes_at_patch = (self.nodes_at_patch > -1).sum(axis=1) any_node_at_patch_closed = (node_status_at_patch == self.BC_NODE_IS_CLOSED).sum( axis=1 ) > (max_nodes_at_patch - 3) absent_patches = any_node_at_patch_closed[self.patches_at_node] bad_patches = np.logical_or(absent_patches, self.patches_at_node == -1) self._patches_present_mask = np.logical_not(bad_patches) self._number_of_patches_present_at_node = np.sum( self._patches_present_mask, axis=1 ) absent_patches = any_node_at_patch_closed[self.patches_at_link] bad_patches = np.logical_or(absent_patches, self.patches_at_link == -1) self._patches_present_link_mask = np.logical_not(bad_patches) self._number_of_patches_present_at_link = np.sum( self._patches_present_link_mask, axis=1 )
[docs] def calc_hillshade_at_node( self, alt=45.0, az=315.0, slp=None, asp=None, unit="degrees", elevs="topographic__elevation", ): """Get array of hillshade. .. codeauthor:: Katy Barnhart <katherine.barnhart@colorado.edu> Parameters ---------- alt : float Sun altitude (from horizon) - defaults to 45 degrees az : float Sun azimuth (CW from north) - defaults to 315 degrees slp : float slope of cells at surface - optional asp : float aspect of cells at surface (from north) - optional (with slp) unit : string 'degrees' (default) or 'radians' - only needed if slp and asp are not provided If slp and asp are both not specified, 'elevs' must be provided as a grid field name (defaults to 'topographic__elevation') or an nnodes-long array of elevation values. In this case, the method will calculate local slopes and aspects internally as part of the hillshade production. Returns ------- ndarray of float Hillshade at each cell. Notes ----- code taken from GeospatialPython.com example from December 14th, 2014 DEJH found what looked like minor sign problems, and adjusted to follow the `ArcGIS algorithm <http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/How_Hillshade_works/009z000000z2000000/>`. Remember when plotting that bright areas have high values. cmap='Greys' will give an apparently inverted color scheme. *cmap='gray'* has white associated with the high values, so is recommended for plotting. Examples -------- >>> import numpy as np >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((5, 5), xy_spacing=1.) >>> z = mg.x_of_node * np.tan(60. * np.pi / 180.) >>> mg.calc_hillshade_at_node(elevs=z, alt=30., az=210.) array([ 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625, 0.625]) LLCATS: NINF SURF """ if slp is not None and asp is not None: if unit == "degrees": (alt, az, slp, asp) = ( np.radians(alt), np.radians(az), np.radians(slp), np.radians(asp), ) elif unit == "radians": if alt > np.pi / 2.0 or az > 2.0 * np.pi: print( "Assuming your solar properties are in degrees, " "but your slopes and aspects are in radians..." ) (alt, az) = (np.radians(alt), np.radians(az)) # ...because it would be super easy to specify radians, # but leave the default params alone... else: raise TypeError("unit must be 'degrees' or 'radians'") elif slp is None and asp is None: if unit == "degrees": (alt, az) = (np.radians(alt), np.radians(az)) elif unit == "radians": pass else: raise TypeError("unit must be 'degrees' or 'radians'") slp, slp_comps = self.calc_slope_at_node(elevs, return_components=True) asp = self.calc_aspect_at_node( slope_component_tuple=slp_comps, unit="radians" ) else: raise TypeError("Either both slp and asp must be set, or neither!") shaded = np.sin(alt) * np.cos(slp) + np.cos(alt) * np.sin(slp) * np.cos( az - asp ) return shaded.clip(0.0)
@property @lru_cache() @make_return_array_immutable def cell_area_at_node(self): """Cell areas in a nnodes-long array. Zeros are entered at all perimeter nodes, which lack cells. Returns ------- ndarray Cell areas as an n_nodes-long array. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((4, 5), xy_spacing=(3, 4)) >>> grid.status_at_node[7] = grid.BC_NODE_IS_CLOSED >>> grid.cell_area_at_node array([ 0., 0., 0., 0., 0., 0., 12., 12., 12., 0., 0., 12., 12., 12., 0., 0., 0., 0., 0., 0.]) LLCATS: CINF NINF CONN """ cell_area_at_node = np.zeros(self.number_of_nodes, dtype=float) cell_area_at_node[self.node_at_cell] = self.area_of_cell return cell_area_at_node
[docs] def reset_status_at_node(self): attrs = [ "_active_link_dirs_at_node", "_status_at_link", "_active_links", "_fixed_links", "_activelink_fromnode", "_activelink_tonode", "_active_faces", "_core_nodes", "_core_cells", "_fixed_links", "_active_adjacent_nodes_at_node", "_fixed_value_boundary_nodes", "_node_at_core_cell", "_link_status_at_node", ] for attr in attrs: try: del self.__dict__[attr] except KeyError: pass try: self.bc_set_code += 1 except AttributeError: self.bc_set_code = 0 try: del self.__dict__["__node_active_inlink_matrix"] except KeyError: pass try: del self.__dict__["__node_active_outlink_matrix"] except KeyError: pass
[docs] def set_nodata_nodes_to_closed(self, node_data, nodata_value): """Make no-data nodes closed boundaries. Sets node status to `BC_NODE_IS_CLOSED` for all nodes whose value of *node_data* is equal to the *nodata_value*. Any links connected to `BC_NODE_IS_CLOSED` nodes are automatically set to `LinkStatus.INACTIVE` boundary. Parameters ---------- node_data : ndarray Data values. nodata_value : float Value that indicates an invalid value. Examples -------- The following example uses the following grid:: *--I--->o------>o------>o ^ ^ ^ ^ I I | | | | | | *--I--->*--I--->o------>o ^ ^ ^ ^ I I I I | | | | *--I--->*--I--->*--I--->* .. note:: Links set to `LinkStatus.ACTIVE` are not shown in this diagram. ``*`` indicates the nodes that are set to `NodeStatus.CLOSED` ``o`` indicates the nodes that are set to `NodeStatus.CORE` ``I`` indicates the links that are set to `LinkStatus.INACTIVE` >>> import numpy as np >>> import landlab as ll >>> mg = ll.RasterModelGrid((3, 4)) >>> mg.status_at_node array([1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1], dtype=uint8) >>> h = np.array([-9999, -9999, -9999, -9999, -9999, -9999, 12345., ... 0., -9999, 0., 0., 0.]) >>> mg.set_nodata_nodes_to_closed(h, -9999) >>> mg.status_at_node array([4, 4, 4, 4, 4, 4, 0, 1, 4, 1, 1, 1], dtype=uint8) LLCATS: BC NINF """ # Find locations where value equals the NODATA code and set these nodes # as inactive boundaries. nodata_locations = np.nonzero(node_data == nodata_value) self.status_at_node[nodata_locations] = self.BC_NODE_IS_CLOSED
[docs] def set_nodata_nodes_to_fixed_gradient(self, node_data, nodata_value): """Make no-data nodes fixed gradient boundaries. Set node status to `BC_NODE_IS_FIXED_VALUE` for all nodes whose value of *node_data* is equal to *nodata_value*. Any links between `BC_NODE_IS_FIXED_GRADIENT` nodes and `BC_NODE_IS_CORE` are automatically set to `LinkStatus.FIXED` boundary status. Parameters ---------- node_data : ndarray Data values. nodata_value : float Value that indicates an invalid value. Examples -------- The following examples use this grid:: *--I--->*--I--->*--I--->*--I--->*--I--->*--I--->*--I--->*--I--->* ^ ^ ^ ^ ^ ^ ^ ^ ^ I I I X X X X X I | | | | | | | | | *--I--->*--I--->*--X--->o o o o o--X--->* ^ ^ ^ ^ ^ ^ ^ ^ ^ I I I | | | | | I | | | | | | | | | *--I--->*--I--->*--X--->o o o o o--X--->* ^ ^ ^ ^ ^ ^ ^ ^ ^ I I I X X X X X I | | | | | | | | | *--I--->*--I--->*--I--->*--I--->*--I--->*--I--->*--I--->*--I--->* .. note:: Links set to `LinkStatus.ACTIVE` are not shown in this diagram. ``X`` indicates the links that are set to `LinkStatus.FIXED` ``I`` indicates the links that are set to `LinkStatus.INACTIVE` ``o`` indicates the nodes that are set to `NodeStatus.CORE` ``*`` indicates the nodes that are set to `BC_NODE_IS_FIXED_GRADIENT` >>> import numpy as np >>> from landlab import RasterModelGrid >>> rmg = RasterModelGrid((4, 9)) >>> rmg.status_at_node # doctest: +NORMALIZE_WHITESPACE array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint8) >>> z = rmg.zeros(at='node') >>> z = np.array([ ... -99., -99., -99., -99., -99., -99., -99., -99., -99., ... -99., -99., -99., 0., 0., 0., 0., 0., -99., ... -99., -99., -99., 0., 0., 0., 0., 0., -99., ... -99., -99., -99., -99., -99., -99., -99., -99., -99.]) >>> rmg.set_nodata_nodes_to_fixed_gradient(z, -99) >>> rmg.status_at_node # doctest: +NORMALIZE_WHITESPACE array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2], dtype=uint8) >>> rmg.status_at_link # doctest: +NORMALIZE_WHITESPACE array([4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 2, 2, 2, 2, 2, 4, 4, 4, 2, 0, 0, 0, 0, 2, 4, 4, 4, 0, 0, 0, 0, 0, 4, 4, 4, 2, 0, 0, 0, 0, 2, 4, 4, 4, 2, 2, 2, 2, 2, 4, 4, 4, 4, 4, 4, 4, 4, 4], dtype=uint8) LLCATS: BC NINF """ # Find locations where value equals the NODATA code and set these nodes # as inactive boundaries. nodata_locations = np.nonzero(node_data == nodata_value) self.status_at_node[nodata_locations] = NodeStatus.FIXED_GRADIENT
@property def unit_vector_sum_xcomponent_at_node(self): """Get array of x-component of unit vector sums at each node. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((3, 3)) >>> len(grid.unit_vector_sum_xcomponent_at_node) == grid.number_of_nodes True >>> grid.unit_vector_sum_xcomponent_at_node array([ 1., 2., 1., 1., 2., 1., 1., 2., 1.]) LLCATS: NINF MEAS """ return self.unit_vector_at_node[:, 0] @property def unit_vector_sum_ycomponent_at_node(self): """Get array of y-component of unit vector sums at each node. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((3, 3)) >>> len(grid.unit_vector_sum_ycomponent_at_node) == grid.number_of_nodes True >>> grid.unit_vector_sum_ycomponent_at_node array([ 1., 1., 1., 2., 2., 2., 1., 1., 1.]) LLCATS: NINF MEAS """ return self.unit_vector_at_node[:, 1]
[docs] def node_is_boundary(self, ids, boundary_flag=None): """Check if nodes are boundary nodes. Check if nodes at given *ids* are boundary nodes. Use the *boundary_flag* to specify a particular boundary type status flag. Parameters ---------- ids : ndarray Node IDs to check. boundary_flag : int, optional A boundary type to check for. Returns ------- ndarray Array of booleans indicating if nodes are boundary nodes. Examples -------- >>> from landlab import RasterModelGrid >>> mg = RasterModelGrid((4, 5)) >>> mg.node_is_boundary([0, 6]) array([ True, False], dtype=bool) >>> mg.node_is_boundary([0, 6], boundary_flag=mg.BC_NODE_IS_CLOSED) array([False, False], dtype=bool) LLCATS: NINF BC """ if boundary_flag is None: return ~(self._node_status[ids] == NodeStatus.CORE) else: return self._node_status[ids] == boundary_flag
[docs] def calc_distances_of_nodes_to_point( self, coord, get_az=None, node_subset=None, out_distance=None, out_azimuth=None ): """Get distances for nodes to a given point. Returns an array of distances for each node to a provided point. If "get_az" is set to 'angles', returns both the distance array and an array of azimuths from up/north. If it is set to 'displacements', it returns the azimuths as a 2xnnodes array of x and y displacements. If it is not set, returns just the distance array. If "node_subset" is set as an ID, or list/array/etc of IDs method returns just the distance (and optionally azimuth) for that node. Point is provided as a tuple (x,y). If out_distance (& out_azimuth) are provided, these arrays are used to store the outputs. This is recommended for memory management reasons if you are working with node subsets. .. note:: Angles are returned in radians but measured clockwise from north. Parameters ---------- coord : tuple of float Coodinates of point as (x, y). get_az: {None, 'angles', 'displacements'}, optional Optionally calculate azimuths as either angles or displacements. The calculated values will be returned along with the distances as the second item of a tuple. node_subset : array_like, optional Calculate distances on a subset of grid nodes. The default is to calculate distances from the provided points to all nodes. out_distance : array_like, optional If provided, put the calculated distances here. Otherwise, create a new array. out_azimuth : array_like, optional If provided, put the calculated distances here. Otherwise, create a new array. Returns ------- ndarray or tuple of ndarray If *get_az* is ``None`` return the array of distances. Otherwise, return a tuple of distances and azimuths. Notes ----- Once you start working with node subsets in Landlab, which can change size between loops, it's quite possible for Python's internal memory management to crap out after large numbers of loops (~>10k). This is to do with the way it block allocates memory for arrays of differing lengths, then cannot free this memory effectively. The solution - as implemented here - is to pre-allocate all arrays as nnodes long, then only work with the first [len_subset] entries by slicing (in a pseudo-C-style). Care has to be taken not to "accidentally" allow Python to allocate a new array you don't have control over. Then, to maintain efficient memory allocation, we create some "dummy" nnode-long arrays to store intermediate parts of the solution in. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((4, 5)) Calculate distances from point at (2., 1.) to a subset of nodes on the grid. >>> grid.calc_distances_of_nodes_to_point((2, 1), ... node_subset=(2, 6, 7, 8, 12)) array([ 1., 1., 0., 1., 1.]) Calculate distances from a point to all nodes on the grid. >>> dist = grid.calc_distances_of_nodes_to_point((2, 1)) >>> dist.shape == (grid.number_of_nodes, ) True >>> dist.take((2, 6, 7, 8, 12)) array([ 1., 1., 0., 1., 1.]) Put the distances into a buffer. >>> out = np.empty(grid.number_of_nodes, dtype=float) >>> dist = grid.calc_distances_of_nodes_to_point((2, 1), ... out_distance=out) >>> out is dist True >>> out.take((2, 6, 7, 8, 12)) array([ 1., 1., 0., 1., 1.]) Calculate azimuths along with distances. The azimuths are calculated in radians but measured clockwise from north. >>> (_, azim) = grid.calc_distances_of_nodes_to_point((2, 1), ... get_az='angles') >>> azim.take((2, 6, 7, 8, 12)) * 180. / np.pi array([ 180., 270., 0., 90., 0.]) >>> (_, azim) = grid.calc_distances_of_nodes_to_point((2, 1), ... get_az='angles', node_subset=(1, 3, 11, 13)) >>> azim * 180. / np.pi array([ 225., 135., 315., 45.]) When calculating displacements, the first row contains displacements in x and the second displacements in y. >>> (_, azim) = grid.calc_distances_of_nodes_to_point((2, 1), ... get_az='displacements', node_subset=(2, 6, 7, 8, 12)) >>> azim array([[ 0., -1., 0., 1., 0.], [-1., 0., 0., 0., 1.]]) LLCATS: NINF MEAS """ if len(coord) != 2: raise ValueError("coordinate must iterable of length 2") if get_az not in (None, "displacements", "angles"): raise ValueError("get_az not understood") if node_subset is not None and np.any(np.isnan(node_subset)): node_subset = None if node_subset is not None: if not isinstance(node_subset, np.ndarray): node_subset = np.array(node_subset) node_subset = node_subset.reshape((-1,)) len_subset = node_subset.size else: len_subset = self.number_of_nodes if out_distance is None: out_distance = np.empty(len_subset, dtype=float) if out_distance.size != len_subset: raise ValueError("output array size mismatch for distances") if get_az is not None: if get_az == "displacements": az_shape = (2, len_subset) else: az_shape = (len_subset,) if out_azimuth is None: out_azimuth = np.empty(az_shape, dtype=float) if out_azimuth.shape != az_shape: raise ValueError("output array mismatch for azimuths") azimuths_as_displacements = np.empty((2, self.number_of_nodes)) dummy_nodes_1 = np.empty(self.number_of_nodes) dummy_nodes_2 = np.empty(self.number_of_nodes) dummy_nodes_3 = np.empty(self.number_of_nodes) if node_subset is None: azimuths_as_displacements[0] = self.node_x - coord[0] azimuths_as_displacements[1] = self.node_y - coord[1] else: azimuths_as_displacements[0, :len_subset] = ( self.node_x[node_subset] - coord[0] ) azimuths_as_displacements[1, :len_subset] = ( self.node_y[node_subset] - coord[1] ) np.square( azimuths_as_displacements[0, :len_subset], out=dummy_nodes_1[:len_subset] ) np.square( azimuths_as_displacements[1, :len_subset], out=dummy_nodes_2[:len_subset] ) np.add( dummy_nodes_1[:len_subset], dummy_nodes_2[:len_subset], out=dummy_nodes_3[:len_subset], ) np.sqrt(dummy_nodes_3[:len_subset], out=out_distance) if get_az: if get_az == "displacements": out_azimuth[:] = azimuths_as_displacements[:, :len_subset] elif get_az == "angles": np.arctan2( azimuths_as_displacements[0, :len_subset], azimuths_as_displacements[1, :len_subset], out=out_azimuth[:len_subset], ) less_than_zero = np.empty(self.number_of_nodes, dtype=bool) np.less(out_azimuth, 0.0, out=less_than_zero[:len_subset]) out_azimuth[less_than_zero[:len_subset]] += 2.0 * np.pi return out_distance, out_azimuth else: return out_distance
@property def all_node_distances_map(self): """Get distances from every node to every other node. Examples -------- >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((3, 4)) >>> distances = grid.all_node_distances_map The shape of the array is ``number_of_nodes`` by ``number_of_nodes`` and distance from a node to itself is zero. >>> distances.shape == (grid.number_of_nodes, grid.number_of_nodes) True >>> distances.diagonal() array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) The distances from the first node to all nodes in its row and all the nodes in its column. >>> distances[0, :4] array([ 0., 1., 2., 3.]) >>> distances[0, ::4] array([ 0., 1., 2.]) LLCATS: NINF MEAS """ if self._all_node_distances_map is None: self._create_all_node_distances_azimuths_maps() return self._all_node_distances_map @property def all_node_azimuths_map(self): """Get azimuths from every node to every other node. Examples -------- >>> import numpy as np >>> from landlab import RasterModelGrid >>> grid = RasterModelGrid((3, 4)) >>> angles = grid.all_node_azimuths_map The shape of the array is ``number_of_nodes`` by ``number_of_nodes`` and azimuth from a node to itself is zero. >>> angles.shape == (grid.number_of_nodes, grid.number_of_nodes) True >>> angles.diagonal() array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) Angles are measured in radians and increase clockwise starting at north. >>> angles *= 180. / np.pi >>> angles[0, :4] array([ 0., 90., 90., 90.]) >>> angles[0, ::4] array([ 0., 0., 0.]) >>> angles[0, ::5] array([ 0., 45., 45.]) LLCATS: NINF MEAS """ if self._all_node_azimuths_map is None: self._create_all_node_distances_azimuths_maps() return self._all_node_azimuths_map def _create_all_node_distances_azimuths_maps(self): """Build distance-azimuth maps. This function creates and stores in the grid field two ``nnodes`` by ``nnodes`` arrays that map the distances and azimuths of all nodes in the grid to all nodes in the grid. This is useful if your module needs to make repeated lookups of distances between the same nodes, but does potentially use up a lot of memory so should be used with caution. The map is symmetrical, so it does not matter whether rows are "from" or "to". The arrays are called: - ``self.all_node_distances_map`` - ``self.all_node_azimuths_map`` Returns ------- tuple of ndarrays Tuple of (distances, azimuths) """ self._all_node_distances_map = np.empty( (self.number_of_nodes, self.number_of_nodes) ) self._all_node_azimuths_map = np.empty( (self.number_of_nodes, self.number_of_nodes) ) node_coords = np.empty((self.number_of_nodes, 2)) node_coords[:, 0] = self.node_x node_coords[:, 1] = self.node_y for i in range(self.number_of_nodes): ( self._all_node_distances_map[i, :], self._all_node_azimuths_map[i, :], ) = self.calc_distances_of_nodes_to_point( (node_coords[i, 0], node_coords[i, 1]), get_az="angles" ) assert np.all(self._all_node_distances_map >= 0.0) return self._all_node_distances_map, self._all_node_azimuths_map # def node_has_boundary_neighbor(self, ids):
[docs] def node_has_boundary_neighbor(self): """Check if ModelGrid nodes have neighbors that are boundary nodes. Checks to see if one of the eight neighbor nodes of node(s) with *id* has a boundary node. Returns True if a node has a boundary node, False if all neighbors are interior. Parameters ---------- ids : int, or iterable of int ID of node to test. Returns ------- boolean ``True`` if node has a neighbor with a boundary ID, ``False`` otherwise. Examples -------- :: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23 >>> from landlab import HexModelGrid >>> grid = HexModelGrid((5, 4)) >>> grid.node_has_boundary_neighbor() array([ True, True, True, True, True, True, True, True, True, True, True, False, False, True, True, True, True, True, True, True, True, True, True, True], dtype=bool) >>> grid.node_has_boundary_neighbor()[6] True >>> grid.node_has_boundary_neighbor()[12] False >>> grid.node_has_boundary_neighbor()[((12, 0),)] array([False, True], dtype=bool) LLCATS: NINF CONN BC """ status_of_neighbor = self._node_status[self.adjacent_nodes_at_node] neighbor_not_core = status_of_neighbor != NodeStatus.CORE bad_neighbor = self.adjacent_nodes_at_node == self.BAD_INDEX neighbor_not_core[bad_neighbor] = False return np.any(neighbor_not_core, axis=1)
# node_has_boundary_neighbor = np.any(neighbor_not_core, axis=1) # return node_has_boundary_neighbor[ids] add_module_functions_to_class(ModelGrid, "mappers.py", pattern="map_*") # add_module_functions_to_class(ModelGrid, 'gradients.py', # pattern='calculate_*') add_module_functions_to_class(ModelGrid, "gradients.py", pattern="calc_*") add_module_functions_to_class(ModelGrid, "divergence.py", pattern="calc_*")