Graph & Connectivity #

Derived connectivity tables, graph traversal, and binary morphology on a mesh.

Everything here is pure: each function reads a Mesh2d and returns new arrays or a new Connectivity, never writing a derived table back onto the mesh it was handed. The edge tables are derived from the faces; when a file declares its own edge_node_connectivity those edges keep their positions (edge-located variables index into them) and any face-implied edge the file omits is appended after them, so a partial edge list still yields a complete face topology.

sparse_adjacency is the shared primitive: connected_components, reverse_cuthill_mckee and the binary morphology all read it.

One morphology pass is simultaneous — the neighbour test is taken from the mask as it stands at the start of the pass, and the border_value imputation joins that same update — so one iteration moves exactly one ring, matching scipy.ndimage.

`pyramids.netcdf.ugrid.graph` #

Mesh graph and connectivity operations for UGRID 2D meshes.

A collection of pure functions that each take a :class:~pyramids.netcdf.ugrid.mesh.Mesh2d (and/or plain numpy arrays) and return plain results — numpy arrays, :class:Connectivity wrappers, scipy sparse matrices, shapely geometries, or a rebuilt Mesh2d for the one topology-transforming helper. None of them mutate their input.

The module groups the graph-flavoured mesh operations behind a single import:

Connectivity builders that pyramids does not derive today: node_node_connectivity, node_face_connectivity, node_edge_connectivity, edge_node_connectivity, edge_face_connectivity, face_edge_connectivity and edge_edge_connectivity.
sparse_adjacency — a scipy sparse adjacency matrix over faces or nodes, the shared primitive the graph algorithms below build on.
Facet remapping (to_node / to_edge / to_face) that averages a per-element field onto another mesh facet via connectivity.
Graph algorithms: connected_components, reverse_cuthill_mckee (+ permute_faces), binary_dilation / binary_erosion, voronoi / circumcenters and the boundary helpers exterior_edges / exterior_faces / bounding_polygon.

scipy (scipy.sparse, scipy.sparse.csgraph, scipy.spatial) and shapely are hard dependencies of the UGRID subsystem, so they are imported at module top.

`node_face_connectivity(mesh)` #

Faces incident to each node (inverse of face_node_connectivity).

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required

Returns:

Name	Type	Description
`Connectivity`	`Connectivity`	Ragged `node -> faces` mapping, `-1` padded.

Source code in src/pyramids/netcdf/ugrid/graph.py

def node_face_connectivity(mesh: Mesh2d) -> Connectivity:
    """Faces incident to each node (inverse of ``face_node_connectivity``).

    Args:
        mesh: The mesh to read.

    Returns:
        Connectivity: Ragged ``node -> faces`` mapping, ``-1`` padded.
    """
    return _invert_connectivity(
        mesh.face_node_connectivity, mesh.n_node, "node_face_connectivity"
    )

`edge_node_connectivity(mesh)` #

The two nodes of each edge, derived from the faces when absent.

A row the file declared but never defined -- one padded with the fill value, or a repeat of an edge already listed -- comes back as a (-1, -1) placeholder. It keeps its position, so an edge-located variable stays aligned with the table, but it describes no edge: :func:_real_edges is the mask of the rows that do, and :func:edge_face_connectivity / :func:edge_edge_connectivity answer an all-fill row at the same positions.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required

Returns:

Name	Type	Description
`Connectivity`	`Connectivity`	`edge -> (node, node)` mapping, with a `(-1, -1)` row wherever the
	`Connectivity`	file's edge list declared an edge it did not define.

Source code in src/pyramids/netcdf/ugrid/graph.py

def edge_node_connectivity(mesh: Mesh2d) -> Connectivity:
    """The two nodes of each edge, derived from the faces when absent.

    A row the file declared but never defined -- one padded with the fill value, or a repeat of
    an edge already listed -- comes back as a ``(-1, -1)`` **placeholder**. It keeps its
    position, so an edge-located variable stays aligned with the table, but it describes no
    edge: :func:`_real_edges` is the mask of the rows that do, and
    :func:`edge_face_connectivity` / :func:`edge_edge_connectivity` answer an all-fill row at
    the same positions.

    Args:
        mesh: The mesh to read.

    Returns:
        Connectivity: ``edge -> (node, node)`` mapping, with a ``(-1, -1)`` row wherever the
        file's edge list declared an edge it did not define.
    """
    edge_node, _, _ = _edge_tables(mesh)
    return Connectivity(
        data=edge_node,
        fill_value=-1,
        cf_role="edge_node_connectivity",
        original_start_index=0,
    )

`node_edge_connectivity(mesh)` #

Edges incident to each node (inverse of edge_node_connectivity).

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required

Returns:

Name	Type	Description
`Connectivity`	`Connectivity`	Ragged `node -> edges` mapping, `-1` padded.

Source code in src/pyramids/netcdf/ugrid/graph.py

def node_edge_connectivity(mesh: Mesh2d) -> Connectivity:
    """Edges incident to each node (inverse of ``edge_node_connectivity``).

    Args:
        mesh: The mesh to read.

    Returns:
        Connectivity: Ragged ``node -> edges`` mapping, ``-1`` padded.
    """
    edge_node, _, _ = _edge_tables(mesh)
    # Keep the edge ids of the real rows only: a placeholder row is not an edge of any node, and
    # its -1 node index would reach np.bincount, which refuses negatives.
    real = _real_edges(edge_node)
    e_ids = np.arange(edge_node.shape[0], dtype=np.intp)[real]
    pairs = edge_node[real]
    keys = np.concatenate([pairs[:, 0], pairs[:, 1]])
    values = np.concatenate([e_ids, e_ids])
    return _group_to_connectivity(keys, values, mesh.n_node, "node_edge_connectivity")

`node_node_connectivity(mesh)` #

Neighbouring nodes of each node, derived from the mesh edges.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required

Returns:

Name	Type	Description
`Connectivity`	`Connectivity`	Ragged `node -> neighbour nodes` mapping, `-1` padded.

Source code in src/pyramids/netcdf/ugrid/graph.py

def node_node_connectivity(mesh: Mesh2d) -> Connectivity:
    """Neighbouring nodes of each node, derived from the mesh edges.

    Args:
        mesh: The mesh to read.

    Returns:
        Connectivity: Ragged ``node -> neighbour nodes`` mapping, ``-1`` padded.
    """
    edge_node, _, _ = _edge_tables(mesh)
    pairs = edge_node[_real_edges(edge_node)]
    keys = np.concatenate([pairs[:, 0], pairs[:, 1]])
    values = np.concatenate([pairs[:, 1], pairs[:, 0]])
    return _group_to_connectivity(keys, values, mesh.n_node, "node_node_connectivity")

`edge_face_connectivity(mesh)` #

The (up to two) faces on either side of each edge.

Edges on the mesh boundary carry a -1 fill on their missing side.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required

Returns:

Name	Type	Description
`Connectivity`	`Connectivity`	`edge -> faces` mapping, `-1` for an exterior side.

Source code in src/pyramids/netcdf/ugrid/graph.py

def edge_face_connectivity(mesh: Mesh2d) -> Connectivity:
    """The (up to two) faces on either side of each edge.

    Edges on the mesh boundary carry a ``-1`` fill on their missing side.

    Args:
        mesh: The mesh to read.

    Returns:
        Connectivity: ``edge -> faces`` mapping, ``-1`` for an exterior side.
    """
    _, edge_face, _ = _edge_tables(mesh)
    return Connectivity(
        data=edge_face,
        fill_value=-1,
        cf_role="edge_face_connectivity",
        original_start_index=0,
    )

`face_edge_connectivity(mesh)` #

The boundary edges of each face.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required

Returns:

Name	Type	Description
`Connectivity`	`Connectivity`	Ragged `face -> edges` mapping, `-1` padded.

Source code in src/pyramids/netcdf/ugrid/graph.py

def face_edge_connectivity(mesh: Mesh2d) -> Connectivity:
    """The boundary edges of each face.

    Args:
        mesh: The mesh to read.

    Returns:
        Connectivity: Ragged ``face -> edges`` mapping, ``-1`` padded.
    """
    _, _, face_edge = _edge_tables(mesh)
    return face_edge

`edge_edge_connectivity(mesh)` #

Edges that share a node with each edge.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required

Returns:

Name	Type	Description
`Connectivity`	`Connectivity`	Ragged `edge -> neighbouring edges` mapping, `-1` padded.

Source code in src/pyramids/netcdf/ugrid/graph.py

def edge_edge_connectivity(mesh: Mesh2d) -> Connectivity:
    """Edges that share a node with each edge.

    Args:
        mesh: The mesh to read.

    Returns:
        Connectivity: Ragged ``edge -> neighbouring edges`` mapping, ``-1`` padded.
    """
    node_edge = node_edge_connectivity(mesh)
    edge_node, _, _ = _edge_tables(mesh)
    n_edge = int(edge_node.shape[0])
    keys_list: list[np.typing.NDArray] = []
    values_list: list[np.typing.NDArray] = []
    for r in range(node_edge.n_elements):
        edges = node_edge.get_element(r)
        if edges.size < 2:
            continue
        left, right = np.meshgrid(edges, edges)
        off_diag = left != right
        keys_list.append(left[off_diag])
        values_list.append(right[off_diag])
    if keys_list:
        keys = np.concatenate(keys_list)
        values = np.concatenate(values_list)
    else:
        keys = np.empty(0, dtype=np.intp)
        values = np.empty(0, dtype=np.intp)
    return _group_to_connectivity(keys, values, n_edge, "edge_edge_connectivity")

`sparse_adjacency(mesh, kind='face')` #

Build a symmetric scipy sparse adjacency matrix for the mesh.

For kind="face" two faces are adjacent when they share an edge (an edge whose edge_face_connectivity has two valid faces). For kind="node" two nodes are adjacent when an edge joins them. The result is the shared primitive consumed by :func:connected_components, :func:reverse_cuthill_mckee and the binary morphology.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required
`kind`	`str`	`"face"` (default) or `"node"`.	`'face'`

Returns:

Type	Description
`Any`	scipy.sparse.csr_matrix: The symmetric `(n, n)` adjacency matrix.

Raises:

Type	Description
`ValueError`	If `kind` is neither `"face"` nor `"node"`.

Source code in src/pyramids/netcdf/ugrid/graph.py

def sparse_adjacency(mesh: Mesh2d, kind: str = "face") -> Any:
    """Build a symmetric scipy sparse adjacency matrix for the mesh.

    For ``kind="face"`` two faces are adjacent when they share an edge (an edge whose
    ``edge_face_connectivity`` has two valid faces). For ``kind="node"`` two nodes are
    adjacent when an edge joins them. The result is the shared primitive consumed by
    :func:`connected_components`, :func:`reverse_cuthill_mckee` and the binary morphology.

    Args:
        mesh: The mesh to read.
        kind: ``"face"`` (default) or ``"node"``.

    Returns:
        scipy.sparse.csr_matrix: The symmetric ``(n, n)`` adjacency matrix.

    Raises:
        ValueError: If ``kind`` is neither ``"face"`` nor ``"node"``.
    """
    if kind == "face":
        _, edge_face, _ = _edge_tables(mesh)
        both_valid = (edge_face[:, 0] != -1) & (edge_face[:, 1] != -1)
        rows = edge_face[both_valid, 0]
        cols = edge_face[both_valid, 1]
        n = mesh.n_face
    elif kind == "node":
        edge_node, _, _ = _edge_tables(mesh)
        pairs = edge_node[_real_edges(edge_node)]
        rows = pairs[:, 0]
        cols = pairs[:, 1]
        n = mesh.n_node
    else:
        raise ValueError(f"kind must be 'face' or 'node', got {kind!r}.")

    sym_rows = np.concatenate([rows, cols])
    sym_cols = np.concatenate([cols, rows])
    data = np.ones(sym_rows.size, dtype=np.int8)
    matrix = sparse.coo_matrix((data, (sym_rows, sym_cols)), shape=(n, n)).tocsr()
    matrix.data[:] = 1
    return matrix

`to_node(mesh, values, source='face')` #

Remap a per-face or per-edge field onto the nodes by averaging.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required
`values`	`NDArray`	Source field, element axis last.	required
`source`	`str`	`"face"` (default) or `"edge"`.	`'face'`

Returns:

Type	Description
`NDArray`	np.ndarray: The node-located field.

Raises:

Type	Description
`ValueError`	If `source` is not `"face"` or `"edge"`.

Source code in src/pyramids/netcdf/ugrid/graph.py

def to_node(
    mesh: Mesh2d, values: np.typing.NDArray, source: str = "face"
) -> np.typing.NDArray:
    """Remap a per-face or per-edge field onto the nodes by averaging.

    Args:
        mesh: The mesh to read.
        values: Source field, element axis last.
        source: ``"face"`` (default) or ``"edge"``.

    Returns:
        np.ndarray: The node-located field.

    Raises:
        ValueError: If ``source`` is not ``"face"`` or ``"edge"``.
    """
    if source == "face":
        conn = node_face_connectivity(mesh)
    elif source == "edge":
        conn = node_edge_connectivity(mesh)
    else:
        raise ValueError(f"source must be 'face' or 'edge', got {source!r}.")
    return _average_over(conn, values)

`to_face(mesh, values, source='node')` #

Remap a per-node or per-edge field onto the faces by averaging.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required
`values`	`NDArray`	Source field, element axis last.	required
`source`	`str`	`"node"` (default) or `"edge"`.	`'node'`

Returns:

Type	Description
`NDArray`	np.ndarray: The face-located field.

Raises:

Type	Description
`ValueError`	If `source` is not `"node"` or `"edge"`.

Source code in src/pyramids/netcdf/ugrid/graph.py

def to_face(
    mesh: Mesh2d, values: np.typing.NDArray, source: str = "node"
) -> np.typing.NDArray:
    """Remap a per-node or per-edge field onto the faces by averaging.

    Args:
        mesh: The mesh to read.
        values: Source field, element axis last.
        source: ``"node"`` (default) or ``"edge"``.

    Returns:
        np.ndarray: The face-located field.

    Raises:
        ValueError: If ``source`` is not ``"node"`` or ``"edge"``.
    """
    if source == "node":
        conn = mesh.face_node_connectivity
    elif source == "edge":
        conn = face_edge_connectivity(mesh)
    else:
        raise ValueError(f"source must be 'node' or 'edge', got {source!r}.")
    return _average_over(conn, values)

`to_edge(mesh, values, source='node')` #

Remap a per-node or per-face field onto the edges by averaging.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required
`values`	`NDArray`	Source field, element axis last.	required
`source`	`str`	`"node"` (default) or `"face"`.	`'node'`

Returns:

Type	Description
`NDArray`	np.ndarray: The edge-located field.

Raises:

Type	Description
`ValueError`	If `source` is not `"node"` or `"face"`.

Source code in src/pyramids/netcdf/ugrid/graph.py

def to_edge(
    mesh: Mesh2d, values: np.typing.NDArray, source: str = "node"
) -> np.typing.NDArray:
    """Remap a per-node or per-face field onto the edges by averaging.

    Args:
        mesh: The mesh to read.
        values: Source field, element axis last.
        source: ``"node"`` (default) or ``"face"``.

    Returns:
        np.ndarray: The edge-located field.

    Raises:
        ValueError: If ``source`` is not ``"node"`` or ``"face"``.
    """
    if source == "node":
        conn = edge_node_connectivity(mesh)
    elif source == "face":
        conn = edge_face_connectivity(mesh)
    else:
        raise ValueError(f"source must be 'node' or 'face', got {source!r}.")
    return _average_over(conn, values)

`connected_components(mesh)` #

Label each face by the connected component it belongs to.

Two faces are in the same component when a path of shared edges joins them — the mesh analogue of a raster cluster labelling.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required

Returns:

Type	Description
`NDArray`	np.ndarray: Per-face component label of length `n_face`.

Source code in src/pyramids/netcdf/ugrid/graph.py

def connected_components(mesh: Mesh2d) -> np.typing.NDArray:
    """Label each face by the connected component it belongs to.

    Two faces are in the same component when a path of shared edges joins them — the mesh
    analogue of a raster cluster labelling.

    Args:
        mesh: The mesh to read.

    Returns:
        np.ndarray: Per-face component label of length ``n_face``.
    """
    adjacency = sparse_adjacency(mesh, kind="face")
    _, labels = csgraph.connected_components(adjacency, directed=False)
    return np.asarray(labels, dtype=np.intp)

`reverse_cuthill_mckee(mesh)` #

Return a face permutation that reduces the face-adjacency bandwidth.

Computes the reverse Cuthill-McKee ordering of the face adjacency. Apply it to the mesh with :func:permute_faces and to every face-located variable with arr[..., perm].

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required

Returns:

Type	Description
`NDArray`	np.ndarray: A permutation of `range(n_face)`.

Source code in src/pyramids/netcdf/ugrid/graph.py

def reverse_cuthill_mckee(mesh: Mesh2d) -> np.typing.NDArray:
    """Return a face permutation that reduces the face-adjacency bandwidth.

    Computes the reverse Cuthill-McKee ordering of the face adjacency. Apply it to the mesh
    with :func:`permute_faces` and to every face-located variable with ``arr[..., perm]``.

    Args:
        mesh: The mesh to read.

    Returns:
        np.ndarray: A permutation of ``range(n_face)``.
    """
    if mesh.n_face == 0:
        # scipy refuses a 0x0 matrix. An empty mesh has an empty ordering, which is the
        # answer rather than an error.
        return np.empty(0, dtype=np.intp)
    adjacency = sparse.csr_matrix(sparse_adjacency(mesh, kind="face"))
    perm = csgraph.reverse_cuthill_mckee(adjacency, symmetric_mode=True)
    return np.asarray(perm, dtype=np.intp)

`permute_faces(mesh, perm)` #

Rebuild a mesh with its faces reordered by perm.

Node coordinates and numbering are unchanged — only the face rows (and any face centre coordinates) are reordered — so the result is congruent to the input. New face k is old face perm[k]; the input mesh is not mutated, and no array is shared with it.

The three face-indexed tables -- face_edge_connectivity, face_face_connectivity and edge_face_connectivity -- are dropped rather than permuted, since their contents index the old face numbering. Derive them again from the result if you need them.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to reorder.	required
`perm`	`NDArray`	A permutation of `range(n_face)`.	required

Returns:

Name	Type	Description
`Mesh2d`	`Mesh2d`	The reordered mesh.

Source code in src/pyramids/netcdf/ugrid/graph.py

def permute_faces(mesh: Mesh2d, perm: np.typing.NDArray) -> Mesh2d:
    """Rebuild a mesh with its faces reordered by ``perm``.

    Node coordinates and numbering are unchanged — only the face rows (and any face centre
    coordinates) are reordered — so the result is congruent to the input. New face ``k`` is
    old face ``perm[k]``; the input mesh is not mutated, and no array is shared with it.

    The three face-indexed tables -- ``face_edge_connectivity``, ``face_face_connectivity``
    and ``edge_face_connectivity`` -- are **dropped** rather than permuted, since their
    contents index the old face numbering. Derive them again from the result if you need
    them.

    Args:
        mesh: The mesh to reorder.
        perm: A permutation of ``range(n_face)``.

    Returns:
        Mesh2d: The reordered mesh.
    """
    perm = np.asarray(perm, dtype=np.intp)
    old_fnc = mesh.face_node_connectivity
    new_fnc = Connectivity(
        data=old_fnc.data[perm].copy(),
        fill_value=old_fnc.fill_value,
        cf_role=old_fnc.cf_role,
        original_start_index=old_fnc.original_start_index,
    )
    face_x = mesh.face_x[perm] if mesh.face_x is not None else None
    face_y = mesh.face_y[perm] if mesh.face_y is not None else None
    return Mesh2d(
        node_x=mesh.node_x.copy(),
        node_y=mesh.node_y.copy(),
        face_node_connectivity=new_fnc,
        # Copy what is carried over, as the permuted arrays above are. Aliasing the
        # untouched tables meant a write through the result reached the input mesh.
        edge_node_connectivity=_copied_connectivity(mesh.edge_node_connectivity),
        face_x=face_x,
        face_y=face_y,
        edge_x=None if mesh.edge_x is None else np.array(mesh.edge_x, copy=True),
        edge_y=None if mesh.edge_y is None else np.array(mesh.edge_y, copy=True),
        crs=mesh.crs,
    )

`binary_dilation(mesh, mask, iterations=1, border_value=False, fixed=None)` #

Grow a True region of a boolean face mask across shared edges.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required
`mask`	`NDArray`	Input boolean face mask (not mutated).	required
`iterations`	`int`	Number of rings to grow.	`1`
`border_value`	`bool`	Value imputed to the mesh exterior (grows inward from the boundary when `True`).	`False`
`fixed`	`NDArray \| None`	Optional boolean faces held at their input value.	`None`

Returns:

Type	Description
`NDArray`	np.ndarray: The dilated boolean face mask.

Source code in src/pyramids/netcdf/ugrid/graph.py

def binary_dilation(
    mesh: Mesh2d,
    mask: np.typing.NDArray,
    iterations: int = 1,
    border_value: bool = False,
    fixed: np.typing.NDArray | None = None,
) -> np.typing.NDArray:
    """Grow a ``True`` region of a boolean face mask across shared edges.

    Args:
        mesh: The mesh to read.
        mask: Input boolean face mask (not mutated).
        iterations: Number of rings to grow.
        border_value: Value imputed to the mesh exterior (grows inward from the boundary
            when ``True``).
        fixed: Optional boolean faces held at their input value.

    Returns:
        np.ndarray: The dilated boolean face mask.
    """
    return _binary_iterate(mesh, mask, True, iterations, border_value, fixed)

`binary_erosion(mesh, mask, iterations=1, border_value=False, fixed=None)` #

Shrink a True region of a boolean face mask across shared edges.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required
`mask`	`NDArray`	Input boolean face mask (not mutated).	required
`iterations`	`int`	Number of rings to shrink.	`1`
`border_value`	`bool`	Value imputed to the mesh exterior (erodes from the boundary when `False`).	`False`
`fixed`	`NDArray \| None`	Optional boolean faces held at their input value.	`None`

Returns:

Type	Description
`NDArray`	np.ndarray: The eroded boolean face mask.

Source code in src/pyramids/netcdf/ugrid/graph.py

def binary_erosion(
    mesh: Mesh2d,
    mask: np.typing.NDArray,
    iterations: int = 1,
    border_value: bool = False,
    fixed: np.typing.NDArray | None = None,
) -> np.typing.NDArray:
    """Shrink a ``True`` region of a boolean face mask across shared edges.

    Args:
        mesh: The mesh to read.
        mask: Input boolean face mask (not mutated).
        iterations: Number of rings to shrink.
        border_value: Value imputed to the mesh exterior (erodes from the boundary when
            ``False``).
        fixed: Optional boolean faces held at their input value.

    Returns:
        np.ndarray: The eroded boolean face mask.
    """
    return _binary_iterate(mesh, mask, False, iterations, border_value, fixed)

`circumcenters(mesh)` #

Per-face circumcenter coordinates for a triangular mesh.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read (every face must be a triangle).	required

Returns:

Type	Description
`NDArray`	np.ndarray: `(n_face, 2)` array of circumcenter `(x, y)` coordinates.

Raises:

Type	Description
`ValueError`	If the mesh is not purely triangular.

Source code in src/pyramids/netcdf/ugrid/graph.py

def circumcenters(mesh: Mesh2d) -> np.typing.NDArray:
    """Per-face circumcenter coordinates for a triangular mesh.

    Args:
        mesh: The mesh to read (every face must be a triangle).

    Returns:
        np.ndarray: ``(n_face, 2)`` array of circumcenter ``(x, y)`` coordinates.

    Raises:
        ValueError: If the mesh is not purely triangular.
    """
    fnc = mesh.face_node_connectivity
    if not fnc.is_triangular():
        raise ValueError("circumcenters is only defined for triangular meshes.")
    tri = fnc.data[:, :3]
    ax, ay = mesh.node_x[tri[:, 0]], mesh.node_y[tri[:, 0]]
    bx, by = mesh.node_x[tri[:, 1]], mesh.node_y[tri[:, 1]]
    cx, cy = mesh.node_x[tri[:, 2]], mesh.node_y[tri[:, 2]]
    d = 2.0 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    a_sq = ax * ax + ay * ay
    b_sq = bx * bx + by * by
    c_sq = cx * cx + cy * cy
    # A degenerate triangle has d == 0 and therefore no circumcentre. Dividing under
    # errstate keeps the inf/nan answer the caller can test for without numpy printing a
    # divide-by-zero or invalid-value warning from inside a library call.
    with np.errstate(divide="ignore", invalid="ignore"):
        ux = (a_sq * (by - cy) + b_sq * (cy - ay) + c_sq * (ay - by)) / d
        uy = (a_sq * (cx - bx) + b_sq * (ax - cx) + c_sq * (bx - ax)) / d
    return np.column_stack([ux, uy])

`voronoi(mesh)` #

Voronoi tessellation of the mesh nodes (the mesh dual).

Uses :class:scipy.spatial.Voronoi over the node coordinates. For a triangular mesh that is Delaunay the Voronoi vertices coincide with the face :func:circumcenters, and interior nodes receive bounded dual cells while boundary nodes receive unbounded (-1) regions.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required

Returns:

Type	Description
`Any`	scipy.spatial.Voronoi: The Voronoi object; `.points` are the node coordinates.

Source code in src/pyramids/netcdf/ugrid/graph.py

def voronoi(mesh: Mesh2d) -> Any:
    """Voronoi tessellation of the mesh nodes (the mesh dual).

    Uses :class:`scipy.spatial.Voronoi` over the node coordinates. For a triangular mesh that
    is Delaunay the Voronoi vertices coincide with the face :func:`circumcenters`, and interior
    nodes receive bounded dual cells while boundary nodes receive unbounded (``-1``) regions.

    Args:
        mesh: The mesh to read.

    Returns:
        scipy.spatial.Voronoi: The Voronoi object; ``.points`` are the node coordinates.
    """
    points = np.column_stack([mesh.node_x, mesh.node_y])
    return Voronoi(points)

`exterior_edges(mesh)` #

Indices of the edges on the mesh boundary.

A boundary edge is one whose edge_face_connectivity has a fill on one side (only one adjacent face).

A dangling edge -- one the file lists that no face implies -- has a fill on both sides, so it is reported here too. That is deliberate: such an edge bounds no face, which makes it part of the mesh's boundary in the only sense this function can measure. A placeholder row (a padded edge the file declared but never defined, or a repeat of an edge already listed) also has two fills, and that one is not reported, because it describes no edge at all.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required

Returns:

Type	Description
`NDArray`	np.ndarray: Boundary edge indices (into the mesh's edge ordering).

Source code in src/pyramids/netcdf/ugrid/graph.py

def exterior_edges(mesh: Mesh2d) -> np.typing.NDArray:
    """Indices of the edges on the mesh boundary.

    A boundary edge is one whose ``edge_face_connectivity`` has a fill on one side (only one
    adjacent face).

    A **dangling** edge -- one the file lists that no face implies -- has a fill on *both*
    sides, so it is reported here too. That is deliberate: such an edge bounds no face, which
    makes it part of the mesh's boundary in the only sense this function can measure. A
    placeholder row (a padded edge the file declared but never defined, or a repeat of an edge
    already listed) also has two fills, and that one is **not** reported, because it describes
    no edge at all.

    Args:
        mesh: The mesh to read.

    Returns:
        np.ndarray: Boundary edge indices (into the mesh's edge ordering).
    """
    edge_node, edge_face, _ = _edge_tables(mesh)
    # A placeholder row (a padded edge the file declared but did not define) touches no
    # face, so it looks exterior on both sides. It is not a boundary edge of the mesh.
    exterior = np.any(edge_face == -1, axis=1) & _real_edges(edge_node)
    return np.flatnonzero(exterior).astype(np.intp)

`exterior_faces(mesh)` #

Indices of the faces that touch the mesh boundary.

A boundary face is one attached to at least one exterior edge.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required

Returns:

Type	Description
`NDArray`	np.ndarray: Sorted unique boundary face indices.

Source code in src/pyramids/netcdf/ugrid/graph.py

def exterior_faces(mesh: Mesh2d) -> np.typing.NDArray:
    """Indices of the faces that touch the mesh boundary.

    A boundary face is one attached to at least one exterior edge.

    Args:
        mesh: The mesh to read.

    Returns:
        np.ndarray: Sorted unique boundary face indices.
    """
    _, edge_face, _ = _edge_tables(mesh)
    exterior = np.any(edge_face == -1, axis=1)
    faces = edge_face[exterior]
    valid = faces[faces != -1]
    return np.unique(valid).astype(np.intp)

`bounding_polygon(mesh)` #

Dissolved outline of the mesh as a shapely (Multi)Polygon, holes preserved.

Unions the face polygons so shared interior edges dissolve away; the result's exterior ring is the mesh outline and any mesh hole appears as an interior ring. Its area equals the mesh footprint. The returned geometry carries no CRS — the caller attaches one.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required

Returns:

Type	Description
`Any`	shapely.Polygon \| shapely.MultiPolygon: The dissolved outline.

Source code in src/pyramids/netcdf/ugrid/graph.py

def bounding_polygon(mesh: Mesh2d) -> Any:
    """Dissolved outline of the mesh as a shapely (Multi)Polygon, holes preserved.

    Unions the face polygons so shared interior edges dissolve away; the result's exterior
    ring is the mesh outline and any mesh hole appears as an interior ring. Its area equals the
    mesh footprint. The returned geometry carries no CRS — the caller attaches one.

    Args:
        mesh: The mesh to read.

    Returns:
        shapely.Polygon | shapely.MultiPolygon: The dissolved outline.
    """
    return unary_union(list(_face_polygons(mesh)))

`value_clusters(mesh, values, lower, upper, nodata=None)` #

Label connected groups of faces whose values fall in [lower, upper] (issue #1379).

The value-aware sibling of :func:connected_components, and the mesh analogue of the raster cluster ("group all the connected values between two bounds"). :func:connected_components labels by topology alone — every face joins a component whatever it holds — so it cannot answer "which pools of flooded cells are there". This restricts the face adjacency to the faces inside the range first, then labels what remains.

Faces outside the range (and no-data faces, which cannot be "in" a range) are excluded from every cluster and labelled 0, so a label is truthy exactly when the face belongs to a cluster. Clusters themselves are numbered 1..k.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required
`values`	`ndarray`	Per-face values, length `n_face`.	required
`lower`	`float`	Inclusive lower bound.	required
`upper`	`float`	Inclusive upper bound.	required
`nodata`	`Any`	Sentinel to exclude, or `None`.	`None`

Returns:

Type	Description
`NDArray`	np.ndarray: Per-face label of length `n_face`; `0` outside every cluster.

Raises:

Type	Description
`ValueError`	`values` is not one per face.

Source code in src/pyramids/netcdf/ugrid/graph.py

def value_clusters(
    mesh: Mesh2d,
    values: np.ndarray,
    lower: float,
    upper: float,
    nodata: Any = None,
) -> np.typing.NDArray:
    """Label connected groups of faces whose values fall in ``[lower, upper]`` (issue #1379).

    The value-aware sibling of :func:`connected_components`, and the mesh analogue of the
    raster ``cluster`` ("group all the connected values between two bounds").
    :func:`connected_components` labels by **topology alone** — every face joins a component
    whatever it holds — so it cannot answer "which pools of flooded cells are there". This
    restricts the face adjacency to the faces inside the range first, then labels what remains.

    Faces outside the range (and no-data faces, which cannot be "in" a range) are excluded from
    every cluster and labelled ``0``, so a label is truthy exactly when the face belongs to a
    cluster. Clusters themselves are numbered ``1..k``.

    Args:
        mesh: The mesh to read.
        values: Per-face values, length ``n_face``.
        lower: Inclusive lower bound.
        upper: Inclusive upper bound.
        nodata: Sentinel to exclude, or ``None``.

    Returns:
        np.ndarray: Per-face label of length ``n_face``; ``0`` outside every cluster.

    Raises:
        ValueError: ``values`` is not one per face.
    """
    array = np.asarray(values)
    if array.ndim != 1 or array.shape[0] != mesh.n_face:
        raise ValueError(
            f"cluster: expected one value per face -- a 1-D array of {mesh.n_face}, got shape "
            f"{array.shape}. Select a single step first, e.g. with isel(time=...)."
        )
    flat = array
    selected = (flat >= lower) & (flat <= upper) & ~nodata_mask(flat, nodata)
    out = np.zeros(mesh.n_face, dtype=np.intp)
    if not selected.any():
        return out
    adjacency = sparse.csr_matrix(sparse_adjacency(mesh, kind="face"))
    # Keep only the edges joining two selected faces, so a cluster cannot leak across a face
    # outside the range. Masking rows and columns by broadcasting avoids `sparse.diags`, whose
    # int8 input would be cast to float64 with a FutureWarning.
    picked = selected.astype(bool)
    restricted = sparse.csr_matrix(
        adjacency.multiply(picked[:, None]).multiply(picked[None, :])
    )
    # `multiply` leaves the masked-out edges in the sparsity structure as explicit zeros, and
    # csgraph treats every *stored* entry as an edge — so without this the excluded faces still
    # joined the clusters they were supposed to separate.
    restricted.eliminate_zeros()
    _, labels = csgraph.connected_components(restricted, directed=False)
    # csgraph gives every isolated node its own component, so the unselected faces each carry a
    # label too. Renumber the selected ones compactly from 1 and leave the rest at 0.
    _, inverse = np.unique(labels[selected], return_inverse=True)
    out[selected] = inverse + 1
    return out

`valid_footprint(mesh, values, nodata=None)` #

The dissolved outline of the faces carrying valid values (issue #1379).

The mesh analogue of the raster footprint ("extract the real coverage of the values"). Distinct from :func:bounding_polygon, which dissolves every face and so describes the mesh domain rather than where data exists; and from polygonize, which dissolves faces of equal value into one polygon per value.

"Valid" means not nodata and not NaN. It deliberately does not consult a variable's active (wet/dry) mask: a dry element has a real value and is still covered by the data, so excluding it would conflate "no data here" with "dry here" — the distinction issue #1331 exists to preserve.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh to read.	required
`values`	`ndarray`	Per-face values, length `n_face`.	required
`nodata`	`Any`	Sentinel to treat as absent, or `None`.	`None`

Returns:

Type	Description
`Any`	shapely.Polygon \| shapely.MultiPolygon: The dissolved coverage, or an empty polygon when
`Any`	no face carries a valid value. Carries no CRS — the caller attaches one.

Raises:

Type	Description
`ValueError`	`values` is not one per face.

Source code in src/pyramids/netcdf/ugrid/graph.py

def valid_footprint(mesh: Mesh2d, values: np.ndarray, nodata: Any = None) -> Any:
    """The dissolved outline of the faces carrying valid values (issue #1379).

    The mesh analogue of the raster ``footprint`` ("extract the real coverage of the values").
    Distinct from :func:`bounding_polygon`, which dissolves **every** face and so describes the
    mesh domain rather than where data exists; and from ``polygonize``, which dissolves faces of
    equal value into one polygon per value.

    "Valid" means not ``nodata`` and not NaN. It deliberately does **not** consult a variable's
    active (wet/dry) mask: a dry element has a real value and is still covered by the data, so
    excluding it would conflate "no data here" with "dry here" — the distinction issue #1331
    exists to preserve.

    Args:
        mesh: The mesh to read.
        values: Per-face values, length ``n_face``.
        nodata: Sentinel to treat as absent, or ``None``.

    Returns:
        shapely.Polygon | shapely.MultiPolygon: The dissolved coverage, or an empty polygon when
        no face carries a valid value. Carries no CRS — the caller attaches one.

    Raises:
        ValueError: ``values`` is not one per face.
    """
    array = np.asarray(values)
    if array.ndim != 1 or array.shape[0] != mesh.n_face:
        raise ValueError(
            f"footprint: expected one value per face -- a 1-D array of {mesh.n_face}, got "
            f"shape {array.shape}. Select a single step first, e.g. with isel(time=...)."
        )
    flat = array
    valid = ~nodata_mask(flat, nodata)
    polygons = [
        polygon
        for index, polygon in enumerate(_face_polygons(mesh))
        if bool(valid[index])
    ]
    if not polygons:
        return shapely.Polygon()
    return unary_union(polygons)

`proximity(mesh, target_faces)` #

Graph distance from every face to the nearest target face.

Uses BFS over the face adjacency graph. A target face has distance zero; an unreachable face (disconnected component with no target) gets -1.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh topology.	required
`target_faces`	`ndarray`	1-D array of face indices that are targets.	required

Returns:

Type	Description
`ndarray`	np.ndarray: Per-face integer distance of length `n_face`.

Examples:

Compute hop distance from face 0 in a two-triangle mesh:

>>> import numpy as np
>>> from pyramids.netcdf.ugrid.connectivity import Connectivity
>>> from pyramids.netcdf.ugrid.mesh import Mesh2d
>>> from pyramids.netcdf.ugrid.graph import proximity
>>> fnc = Connectivity(
...     np.array([[0, 1, 2], [1, 3, 2]]),
...     -1, "face_node_connectivity", 0,
... )
>>> mesh = Mesh2d(
...     node_x=np.array([0.0, 1.0, 0.5, 1.5]),
...     node_y=np.array([0.0, 0.0, 1.0, 1.0]),
...     face_node_connectivity=fnc,
... )
>>> proximity(mesh, np.array([0]))
array([0, 1])

All faces as targets gives distance zero everywhere:

>>> import numpy as np
>>> from pyramids.netcdf.ugrid.connectivity import Connectivity
>>> from pyramids.netcdf.ugrid.mesh import Mesh2d
>>> from pyramids.netcdf.ugrid.graph import proximity
>>> fnc = Connectivity(
...     np.array([[0, 1, 2], [1, 3, 2]]),
...     -1, "face_node_connectivity", 0,
... )
>>> mesh = Mesh2d(
...     node_x=np.array([0.0, 1.0, 0.5, 1.5]),
...     node_y=np.array([0.0, 0.0, 1.0, 1.0]),
...     face_node_connectivity=fnc,
... )
>>> proximity(mesh, np.array([0, 1]))
array([0, 0])

`sieve(mesh, values, threshold)` #

Remove small connected regions below a size threshold.

Faces in a connected component smaller than threshold are set to NaN. The result is always float64 regardless of the input dtype.

Parameters:

Name	Type	Description	Default
`mesh`	`Mesh2d`	The mesh topology.	required
`values`	`ndarray`	1-D variable data (one value per face).	required
`threshold`	`int`	Minimum component size to keep.	required

Returns:

Type	Description
`ndarray`	np.ndarray: A copy of values with small-region faces blanked.

Examples:

Blank faces in a two-face mesh (one component of size 2) when the threshold exceeds the component size:

>>> import numpy as np
>>> from pyramids.netcdf.ugrid.connectivity import Connectivity
>>> from pyramids.netcdf.ugrid.mesh import Mesh2d
>>> from pyramids.netcdf.ugrid.graph import sieve
>>> fnc = Connectivity(
...     np.array([[0, 1, 2], [1, 3, 2]]),
...     -1, "face_node_connectivity", 0,
... )
>>> mesh = Mesh2d(
...     node_x=np.array([0.0, 1.0, 0.5, 1.5]),
...     node_y=np.array([0.0, 0.0, 1.0, 1.0]),
...     face_node_connectivity=fnc,
... )
>>> sieve(mesh, np.array([5.0, 7.0]), threshold=3)
array([nan, nan])

Keep both faces when the threshold is small enough:

>>> import numpy as np
>>> from pyramids.netcdf.ugrid.connectivity import Connectivity
>>> from pyramids.netcdf.ugrid.mesh import Mesh2d
>>> from pyramids.netcdf.ugrid.graph import sieve
>>> fnc = Connectivity(
...     np.array([[0, 1, 2], [1, 3, 2]]),
...     -1, "face_node_connectivity", 0,
... )
>>> mesh = Mesh2d(
...     node_x=np.array([0.0, 1.0, 0.5, 1.5]),
...     node_y=np.array([0.0, 0.0, 1.0, 1.0]),
...     face_node_connectivity=fnc,
... )
>>> sieve(mesh, np.array([5.0, 7.0]), threshold=1)
array([5., 7.])