Meshing#

Mesh#

Tools for turning data into meshes, then writing them to file.

Works with and assumes all polyhedra cells (and polygon faces).

Requires building a reasonably recent version of Exodus to get the associated exodus.py wrappers in python3.

Note that this is typically done in your standard ATS installation, assuming you have built your Amanzi TPLs with shared libraries (the default through bootstrap).

In that case, simply ensure that ${AMANZI_TPLS_DIR}/SEACAS/lib is in your PYTHONPATH.

class watershed_workflow.mesh.Edge(*args)[source]#

A pair of vertex indices forming a mesh edge.

Edges are order-independent: Edge(i, j) == Edge(j, i).

Parameters:

*args (int or Tuple[int, int]) – Either two ints Edge(i, j) or a single tuple Edge((i, j)).

class watershed_workflow.mesh.LabeledSet(name: str, setid: int, entity: str, ent_ids: List[Any], to_extrude: bool = False)[source]#

A generic collection of entities.

validate(size: int, is_tuple: bool = False) → None[source]#

Validate the labeled set.

Parameters:

  • size (int) – Expected size for validation.

  • is_tuple (bool, optional) – Whether entities are tuples (edges). Default is False.

Raises:

AssertionError – If validation fails.

class watershed_workflow.mesh.Mesh2D(coords, conn, labeled_sets: List[LabeledSet] = NOTHING, crs: CRS | None = None, cell_data: DataFrame | None = None, eps: float = 0.001, check_handedness: bool = True, validate: bool = False)[source]#

A 2D mesh class.

Parameters:

  • coords (np.ndarray(NCOORDS,NDIMS)) – Array of coordinates of the 2D mesh.

  • conn (list(lists)) – List of lists of indices into coords that form the cells.

  • labeled_sets (list(LabeledSet), optional) – List of labeled sets to add to the mesh.

  • crs (CRS) – Keep this as a property for future reference.

  • cell_data (Optional[pandas.DataFrame]) – A DataFrame of length len(conn) including assorted cell-based data.

  • eps (float, optional=0.01) – A small measure of length between coords.

  • check_handedness (bool, optional=True) – If true, ensure all cells are oriented so that right-hand rule ordering of the vertices points up.

  • validate (bool, optional=False) – If true, validate coordinates and connections post-construction.

  • Note ((coords, conn) may be output provided by a)

  • call. (watershed_workflow.triangulation.triangulate())

property boundary_edges#

Return edges in the boundary of the mesh, ordered around the boundary.

property cell_edges#

A map from cell to list of edges.

property cell_to_cells#

A list of length ncells, each entry is a list of neighboring cells.

property centroids#

Cell centroids.

checkHandedness(conn=None) → List[int] | None[source]#

Ensures all cells are oriented via the right-hand-rule, i.e. in the +z direction.

clearGeometryCache() → None[source]#

If coordinates are changed, any computed, cached geometry must be recomputed. It is the USER’s responsibility to call this function if any coords are changed!

computeCentroid(c) → ndarray[source]#

Computes, based on coords, the centroid of a cell with ID c.

Note this ALWAYS recomputes the value, not using the cache.

property conn#

Note that conn is immutable because changing this breaks all the other properties, which may be cached. To change topology, one must construct a new mesh!

property dim: int#

Spatial dimension of the mesh.

Returns:

Number of spatial dimensions.

Return type:

int

property edge_cells#

A map from edge to lists of cells that edge touches.

property edge_centroids#

Edge centroids.

classmethod from_Transect(x, z, width=1, **kwargs)[source]#

Creates a 2D surface strip mesh from transect data

getNextAvailableLabeledSetID() → int[source]#

Returns next available LS id.

property num_cells: int#

Number of cells in the mesh.

Returns:

Number of cells.

Return type:

int

property num_edges: int#

Number of edges in the mesh.

Returns:

Number of edges.

Return type:

int

property num_vertices: int#

Number of vertices in the mesh.

Returns:

Number of vertices.

Return type:

int

partition(nparts: int, reorder: bool = True) → Mesh2D[source]#

Partitions the mesh, adding a cell_data column for partition number.

Parameters:

  • nparts (int) – Number of parts to partition cells into.

  • reorder (bool, optional) – If True, also creates a new mesh in the partition ordering. Default is True.

Returns:

If reorder, the partitioned and reordered mesh; otherwise returns self, which includes partition info in cell_data.

Return type:

Mesh2D

plot(facecolors=None, ax=None, cmap=None, vmin=None, vmax=None, norm=None, colorbar=True, **kwargs) → PolyCollection[source]#

Plot the flattened 2D mesh.

plotVertices(vertex_values, ax=None, cmap=None, vmin=None, vmax=None, norm=None, s=None, colorbar=True, label=None, **kwargs)[source]#

Plot vertex-based data as a scatter plot.

Parameters:

  • vertex_values (array-like) – Values at each vertex (length must equal num_vertices)

  • ax (matplotlib.Axes, optional) – Axes to plot on. If None, creates new figure/axes.

  • cmap (str or matplotlib.colors.Colormap, optional) – Colormap for mapping values to colors. Default is ‘viridis’.

  • vmin (float, optional) – Min and max values for color mapping. If None, uses data range.

  • vmax (float, optional) – Min and max values for color mapping. If None, uses data range.

  • norm (matplotlib.colors.Normalize, optional) – Normalization for color mapping. If None, uses linear normalization.

  • s (float, optional) – Marker size. If None, automatically computed based on mesh size.

  • colorbar (bool, optional) – If True, adds colorbar. Default is True.

  • label (str, optional) – Label for colorbar.

  • **kwargs (dict) – Additional keyword arguments passed to scatter()

Returns:

The scatter plot collection

Return type:

PathCollection

classmethod read_VTK(filename)[source]#

Constructor from a VTK file.

classmethod read_VTK_Simplices(filename)[source]#

Constructor from an structured VTK file.

classmethod read_VTK_Unstructured(filename)[source]#

Constructor from an unstructured VTK file.

reorder(new_cell_order) → Mesh2D[source]#

Creates a new mesh with reordered cells.

Parameters:

new_order (list[int]) – A list of length num_cells, indicating the new ordering of existing cell indices. The new mesh’s cell 0 will be self’s new_cell_order[0] cell.

Returns:

  • Mesh2D – The reordered mesh.

  • Note this also tries to reorder the vertex coordinates in a

  • sane way, by reiterating over cell vertices and renumbering.

  • Edges are naturally reordered in the new mesh.

to_dataframe(include_labeled_sets: bool = False) → GeoDataFrame[source]#

Convert the mesh to a GeoDataFrame with each cell as a row.

to_dual()[source]#

Creates a 2D surface mesh from a primal 2D triangular surface mesh.

Returns:

  • dual_vertices (np.ndarray)

  • dual_conn (list(lists)) – The vertices and cell_to_vertex_conn of a 2D, polygonal, nearly Voronoi mesh that is the truncated dual of self. Here we say nearly because the boundary triangles (see note below) are not Voronoi.

  • dual_from_primal_mapping (np.array( (len(dual_conn),), ‘i’)) – Mapping from a dual cell to the primal vertex it is based on. Without the boundary, this would be simply arange(0,len(dual_conn)), but beacuse we added triangles on the boundary, this mapping is useful.

Note

  • Vertices of the dual are on circumcenters of the primal triangles.

  • Centers of the dual are numbered by vertices of the primal mesh, modulo the boundary.

  • At the boundary, the truncated dual polygon may be non-convex. To avoid this issue, we triangulate boundary polygons.

We do not simply construct a mesh because likely the user needs to elevate, generate material IDs, labeled sets, etc, prior to mesh construction.

transform(mat: ndarray | None = None, shift: ndarray | None = None) → None[source]#

Transform a 2D mesh.

Parameters:

  • mat (np.ndarray, optional) – 3x3 transformation matrix. Default is identity.

  • shift (np.ndarray, optional) – 3-element translation vector. Default is zero.

writeVTK(filename: str) → None[source]#

Writes to VTK.

Parameters:

filename (str) – Output VTK filename.

class watershed_workflow.mesh.Mesh3D(coords, face_to_vertex_conn, cell_to_face_conn, labeled_sets: List[LabeledSet] = NOTHING, side_sets: List[SideSet] = NOTHING, material_ids: ndarray = None, crs: CRS | None = None, cell_data: DataFrame | None = None, eps: float = 0.001, validate: bool = False)[source]#

A 3D mesh class.

Parameters:

  • coords (np.ndarray(NCOORDS,3)) – Array of coordinates of the 3D mesh.

  • face_to_vertex_conn (list(lists)) – List of lists of indices into coords that form the faces.

  • cell_to_face_conn (list(lists)) – List of lists of indices into face_to_vertex_conn that form the cells.

  • side_sets (list(SideSet), optional) – List of side sets to add to the mesh.

  • labeled_sets (list(LabeledSet), optional) – List of labeled sets to add to the mesh.

  • material_ids (np.array((len(cell_to_face_conn),),'i'), optional) – Array of length num_cells that specifies material IDs

  • crs (CRS) – Keep the coordinate system for reference.

  • cell_data (dict | pandas.DataFrame) – Extra cell-based data, stored as a dataframe.

  • eps (float, optional=0.01) – A small measure of length between coords.

  • (coords (Note that)

  • a (conn) may be output provided by)

  • call. (watershed_workflow.triangulation.triangulate())

classmethod extruded_Mesh2D(mesh2D, layer_types, layer_data, ncells_per_layer, mat_ids)[source]#

Uniformly extrude a 2D mesh to make a 3D mesh.

Layers of potentially multiple sets of cells are extruded downward in the vertical. The cell dz is uniform horizontally and vertically through the layer in all but the ‘snapped’ case, where it is uniform vertically but not horizontally.

Parameters:

  • mesh2D (Mesh2D) – The 2D mesh to extrude

  • layer_types (str, list[str]) – One of [‘snapped’, ‘function’, ‘vertex’, ‘cell’,] or a list of these strings, describing the extrusion method for each layer. If only a single string is supplied, all layers are of this type. See layer_data below.

  • layer_data (list[data]) –

    Data required, one for each layer. The data is specific to the layer type:

    • ’constant’ : float, layer thickness

    • ’snapped’ : float, the bottom z-coordinate of the layer. Cells are extruded uniformly to whatever thickness required to match this bottom coordinate.

    • ’function’ : function, layer_thickness = func(x,y)

    • ’vertex’ : np.array((mesh2D.num_vertices,), float) layer thickness for each vertex

    • ’cell’ : np.array((mesh2D.num_cells,), float) interpolates nodal layer thickness from neighboring cell thicknesses

  • ncells_per_layer (int, list[int]) – Either a single integer (same number of cells in all layers) or a list of number of cells in each layer

  • mat_ids (int, list[int], np.ndarray((num_layers, mesh2D.num_cells), dtype=int)) – Material ID for each cell in each layer. If an int, all cells in all layers share thesame material ID. If a list or 1D array, one material ID per layer. If a 2D array, provides all material IDs explicitly.

Returns:

The extruded, 3D mesh.

Return type:

Mesh3D

getNextAvailableLabeledSetID() → int[source]#

Returns next available LS id.

static summarizeExtrusion(layer_types, layer_data, ncells_per_layer, mat_ids, surface_cell_id=0) → None[source]#

Summarizes extruded data by printing info to log file.

This is useful in rapidly debugging and understanding the layering before you do the extrusion process.

validate() → None[source]#

Checks the validity of the mesh, or throws an AssertionError.

writeExodus(filename: str, element_block_mode: str) → None[source]#

Write the 3D mesh to ExodusII using arbitrary polyhedra spec.

Parameters:

  • filename (str) – Output Exodus filename.

  • element_block_mode (str, optional) –

    • “one block” (default): preserves sequential columnar ordering of elements by writing all elements to a single element block, and treating material IDs as element sets.

    • ”material id” (old method): one element block per material ID – reorders cells by material ID.

writeVTK(filename: str) → None[source]#

Writes to VTK.

Note, this just writes the topology/geometry information, for WEDGE type meshes (extruded triangles). No labeled sets are written. Prefer to use writeExodus() for a fully featured mesh.

Parameters:

filename (str) – Output VTK filename.

class watershed_workflow.mesh.SideSet(name: str, setid: int, cell_list: list[int], side_list: list[int])[source]#

A collection of faces in cells.

validate(cell_faces: List[List[int]]) → None[source]#

Validate the side set against cell faces.

Parameters:

cell_faces (List[List[int]]) – List of cell face connections.

Raises:

AssertionError – If validation fails.

watershed_workflow.mesh.cache(func: Callable[[Any], Any]) → Callable[[Any], Any][source]#

A caching decorator for instance methods.

This decorator caches the result of a method call in an instance attribute. Note: Only works with __dict__ classes, not slotted classes.

Parameters:

func (Callable[[Any], Any]) – The function to cache.

Returns:

The wrapped function with caching.

Return type:

Callable[[Any], Any]

watershed_workflow.mesh.computeTelescopeFactor(ncells: int, dz: float, layer_dz: float) → float[source]#

Calculates a telescoping factor to fill a given layer.

Calculates a constant geometric factor, such that a layer of thickness layer_dz is perfectly filled by ncells in the vertical, where the top cell is dz in thickness and each successive cell grows by a factor of that factor.

Parameters:

  • ncells (int) – Number of cells (in the vertical) needed.

  • dz (float) – Top cell’s thickness in the vertical.

  • layer_dz (float) – Thickness of the total layer.

Returns:

The telescoping factor.

Return type:

float

watershed_workflow.mesh.createSubmesh(m2: Mesh2D, shp: BaseGeometry) → Tuple[Dict[int, int], Dict[int, int], Mesh2D][source]#

Given a shape that contains some cells of m2, create the submesh.

watershed_workflow.mesh.mergeMeshes(meshes: List[Mesh2D]) → Mesh2D[source]#

Combines multiple 2D meshes into a single mesh.

It is assumed that the meshes to be combined have common vertices on the shared edge (no steiner points). labeledsets should be added after merging the mesh. Option of merging pre-existing labeledsets in the to-be-merged meshes will be added soon

Parameters:

meshes (list(mesh.Mesh2D)) – The list of meshes to be merged

Returns:

combined mesh

Return type:

mesh.Mesh2D

watershed_workflow.mesh.mergeTwoMeshes(mesh1: Mesh2D, mesh2: Mesh2D) → Mesh2D[source]#

merge two meshes (mesh.Mesh2D objects)

watershed_workflow.mesh.optimizeDzs(dz_begin: float, dz_end: float, thickness: float, num_cells: int, p_thickness: float = 1000, p_dz: float = 10000, p_increasing: float = 1000, p_smooth: float = 10, tol: float = 1) → Tuple[ndarray, float][source]#

Tries to optimize dzs

watershed_workflow.mesh.refineCorridorTriangles(m2: Mesh2D, river_corrs: List[Polygon]) → Mesh2D[source]#

Given a mesh, refine all triangles where all three vertices are on the river corridor.

This deals with both interior junction triangles and sharp-angle reach triangles.

watershed_workflow.mesh.refineTriangle(m2: Mesh2D, c: int) → Mesh2D[source]#

Make a new mesh by refining a triangular cell.

Note that cell c must be: - a triangle - that is not on the boundary - whose neighboring cells are also triangles

watershed_workflow.mesh.refineTriangles(m2: Mesh2D, to_refine: List[int]) → Mesh2D[source]#

Refine a set of triangles, making a new mesh.

Note that cell in to_refine must be:

  • a triangle

  • that is not on the boundary

  • whose neighboring cells are also triangles

  • distinct – for each pair of cells c1,c2 in to_refine, the set of c1 and its neighbors must not overlap with the set of c2 and its neighbors.

River-aligned Mesh#

creates river mesh using quad, pentagon and hexagon elements

watershed_workflow.river_mesh.adjustHUCsToRiverMesh(hucs: SplitHUCs, river: River, coords: ndarray) → None[source]#

Adjust HUC segments that touch reach endpoints to match the corridor coordinates.

Parameters:

  • hucs (SplitHUCs) – Split HUCs object to adjust.

  • river (River) – River network with mesh coordinates.

  • coords (np.ndarray) – Array of mesh coordinates.

watershed_workflow.river_mesh.computeLine(p1: ndarray, p2: ndarray) → Tuple[float, float, float][source]#

Compute line coefficients (Ax + By + C = 0) for a line defined by two points.

Parameters:

  • p1 (np.ndarray) – First point coordinates.

  • p2 (np.ndarray) – Second point coordinates.

Returns:

Line coefficients (A, B, C) such that Ax + By + C = 0.

Return type:

Tuple[float, float, float]

watershed_workflow.river_mesh.createRiverMesh(river: River, computeWidth: Callable[[River], float], elems_gid_start: int = 0, check_convexity: bool = True)[source]#

Returns list of elems and river corridor polygons for a given list of river trees

Parameters:

  • river (River) – River tree along which river mesh is to be created.

  • computeWidth (Callable[[River, ], float]) – Function that computes the width for a given reach.

  • elems_gid_start (int, optional) – Starting global ID for elements, by default 0.

  • check_convexity (bool, optional) – If True, check each element for convexity and attempt to fix non-convex tip elements after projection. Set to False to skip this pass entirely, which is useful when deliberately testing degenerate width configurations that would otherwise raise inside the convexity fixer. Default is True.

Returns:

  • corrs (list(shapely.geometry.Polygon)) – List of river corridor polygons, one per river, storing the coordinates used in elems.

  • elems (list(list)) – List of river elements, each element a list of indices into corr.coords.

watershed_workflow.river_mesh.createRiversMesh(hucs: SplitHUCs, rivers: List[River], computeWidth: Callable[[River], float], ax: Axes | None = None, plot: bool = False) → Tuple[ndarray, List[List[int]], List[Polygon], List[Point], GeoDataFrame | None][source]#

Create meshes for each river and merge them.

Parameters:

  • hucs (SplitHUCs) – Split HUCs object for mesh adjustment.

  • rivers (List[River]) – List of river networks to mesh.

  • computeWidth (Callable[[River], float]) – Function to compute the river width for each reach (given as a River object). This callable can either return a constant value, or dynamically fetch a value based on stream order, properties, or a user-defined rule.

  • ax (matplotlib.axes.Axes, optional) – Axes for debugging plots, by default None.

Returns:

Tuple containing coordinates, elements, corridors, hole points, and intersections dataframe.

Return type:

Tuple[np.ndarray, List[List[int]], List[shapely.geometry.Polygon], List[shapely.geometry.Point], gpd.GeoDataFrame | None]

watershed_workflow.river_mesh.findIntersection(line1: Tuple[float, float, float], line2: Tuple[float, float, float], debug: bool = False) → ndarray | None[source]#

Find the intersection point of two lines given by coefficients.

Parameters:

  • line1 (Tuple[float, float, float]) – First line coefficients (A, B, C).

  • line2 (Tuple[float, float, float]) – Second line coefficients (A, B, C).

  • debug (bool, optional) – Whether to print debug information, by default False.

Returns:

Intersection point coordinates, or None if lines are parallel.

Return type:

np.ndarray | None

watershed_workflow.river_mesh.fixConvexity(reach: River, e_coords: ndarray, computeWidth: Callable[[River], float]) → ndarray[source]#

Snap element coordinates onto the convex hull while respecting upstream stream width.

Parameters:

  • reach (River) – River reach containing the element.

  • e_coords (np.ndarray) – Element coordinates to fix.

  • computeWidth (Callable[[River], float]) – Function to compute river width.

Returns:

Fixed element coordinates.

Return type:

np.ndarray

watershed_workflow.river_mesh.projectJunction(reach: River, child_idx: int, computeWidth: Callable[[River], float], debug: bool = False) → ndarray[source]#

Find points around the junction between reach and its children.

Parameters:

  • reach (River) – Parent reach containing the junction.

  • child_idx (int) – Index of the child reach at the junction.

  • computeWidth (Callable[[River], float]) – Function to compute river width.

  • debug (bool, optional) – Whether to print debug information, by default False.

Returns:

Junction point coordinates.

Return type:

np.ndarray

watershed_workflow.river_mesh.projectOne(p_up: ndarray, p: ndarray, width: float) → ndarray[source]#

Find a point p_out that is width away from p and such that p_up –> p is right-perpendicular to p –> p_out.

Parameters:

  • p_up (np.ndarray) – Upstream point coordinates.

  • p (np.ndarray) – Reference point coordinates.

  • width (float) – Distance to project perpendicular to the line.

Returns:

Projected point coordinates.

Return type:

np.ndarray

watershed_workflow.river_mesh.projectTwoBisector(p_up: ndarray, p: ndarray, p_dn: ndarray, width: float, debug: bool = False) → ndarray[source]#

Find a point that is width away from p along the bisector of p_up –> p –> p_dn.

Unlike projectTwoMiter(), which uses a miter approach and accepts two widths, this uses the angle bisector of the two segments and a single width. The result is always exactly width away from p, so it cannot overshoot regardless of the bend angle or width ratio between adjacent reaches.

When the two segments are antiparallel (straight line, zero bend), the bisector is degenerate and the function falls back to a simple perpendicular projection via projectOne() using the upstream segment.

Parameters:

  • p_up (np.ndarray) – Upstream point coordinates.

  • p (np.ndarray) – Middle point coordinates.

  • p_dn (np.ndarray) – Downstream point coordinates.

  • width (float) – Distance from p to the projected bank point, measured along the bisector.

  • debug (bool, optional) – Whether to log debug information, by default False.

Returns:

Projected bank point coordinates.

Return type:

np.ndarray

watershed_workflow.river_mesh.projectTwoClampedMiter(p_up: ndarray, p: ndarray, p_dn: ndarray, width1: float, width2: float, debug: bool = False) → ndarray[source]#

Find a bank point using the miter approach, clamped to avoid overshoot.

Computes the miter intersection via _projectTwoMiter(), then applies one of two strategies depending on the bend angle theta at p:

  • theta > 90 deg (nearly-straight reach, p_up-p-p_dn angle close to 180): the miter direction becomes unreliable; fall back to the bisector direction with the average width, similarly capped.

  • theta <= 90 deg (sharp bend): the miter direction is reliable; cap |m - p| to 3/4 * min(|p_up - p|, |p_dn - p|) to prevent overshoot.

In both cases the fallback also triggers when |m - p| already exceeds the cap (e.g. parallel segments where _projectTwoMiter returns a perpendicular).

Parameters:

  • p_up (np.ndarray) – Upstream point coordinates.

  • p (np.ndarray) – Middle point coordinates.

  • p_dn (np.ndarray) – Downstream point coordinates.

  • width1 (float) – Width for upstream segment.

  • width2 (float) – Width for downstream segment.

  • debug (bool, optional) – Whether to log debug information, by default False.

Returns:

Projected bank point coordinates.

Return type:

np.ndarray

watershed_workflow.river_mesh.translateLinePerpendicular(line: Tuple[float, float, float], distance: float) → Tuple[float, float, float][source]#

Translate a line by a specified distance in the direction perpendicular to the line.

Parameters:

  • line (Tuple[float, float, float]) – Tuple of line coefficients (A, B, C).

  • distance (float) – Scalar distance to translate the line.

Returns:

Tuple of new line coefficients (A, B, C).

Return type:

Tuple[float, float, float]

Triangulation#

Triangulates polygons

class watershed_workflow.triangulation.Nodes(decimals: int = 3)[source]#

A collection of nodes that are indexed in the order they are added.

Note this uses round() for an efficient solution, which is potentially fragile if numbers were not generated consistently. In this use case, however, it should be safe – numbers were originally rounded (in watershed_workflow.config), and then any function which inserted new points or moved points were always careful to ensure that all duplicates were always assigned the identical values – i.e. all math was done BEFORE assignment. So duplicates have identical representations in floating point. This suggests we shouldn’t have to round here at all, but we do anyway to keep things cleaner.

class watershed_workflow.triangulation.NodesEdges(objlist: List[LineString | Polygon] | None = None)[source]#

A collection of nodes and edges.

add(obj: LineString | Polygon) → None[source]#

Adds nodes and edges from obj into collection.

check(tol: float) → None[source]#

Checks consistency of the internal representation.

watershed_workflow.triangulation.connectOnewayTrip(inds: List[int]) → List[Tuple[int, int]][source]#

Connect indices in edges in a oneway fashion

watershed_workflow.triangulation.connectRoundTrip(inds: List[int]) → List[Tuple[int, int]][source]#

Connect indices in edges in a round-trip fashion

watershed_workflow.triangulation.orient(e: Tuple[int, int]) → Tuple[int, int] | None[source]#

Orient an edge consistently.

Parameters:

e (Tuple[int, int]) – Edge as tuple of two vertex indices.

Returns:

Oriented edge or None if vertices are identical.

Return type:

Optional[Tuple[int, int]]

watershed_workflow.triangulation.refineByMaxArea(max_area)[source]#

Returns a refinement function based on max area, for use with Triangle.

watershed_workflow.triangulation.refineByMaxEdgeLength(edge_length)[source]#

Returns a refinement function based on max edge length, for use with Triangle.

watershed_workflow.triangulation.refineByPolygons(polygons, areas)[source]#

Returns a graded refinement function based upon polygon area limits, for use with Triangle.

Triangle area must be smaller than the area limit for the polygon when the triangle centroid is within the polygon.

watershed_workflow.triangulation.refineByRiverDistance(near_distance, near_area, away_distance, away_area, rivers)[source]#

Returns a graded refinement function based upon a distance function from rivers, for use with Triangle.

Triangle area must be smaller than near_area when the triangle centroid is within near_distance from the river network. Area must be smaller than away_area when the triangle centroid is at least away_distance from the river network. Area must be smaller than a linear interpolant between near_area and away_area when between near_distance and away_distance from the river network.

watershed_workflow.triangulation.refineByStreamTriangles(river_corrs)[source]#

Returns a refinement function for triangles that have all three vertices on stream mesh.

watershed_workflow.triangulation.triangulate(hucs: SplitHUCs, internal_boundaries: List[LineString] | None = None, hole_points: List[Point] | None = None, additional_vertices: List[Tuple[float, float]] | None = None, tol: float = 1.0, **kwargs) → Tuple[ndarray, ndarray][source]#

Triangulates HUCs and rivers.

Note, refinement of a given triangle is done if any of the provided criteria is met.

Parameters:

  • hucs (SplitHUCs) – A split-form HUC object from, e.g., get_split_form_hucs()

  • internal_boundaries (list, optional) – List of shapely objects or RiverTrees or other iterable collections of coordinates used as internal boundaries that must be included in the mesh.

  • hole_points (list(shapely.Point), optional) – List of points inside the polygons to be left as holes/voids (excluded from mesh).

  • additional_vertices (list(Tuple[float, float]), optional) – List of points to be included in the triangulation.

  • tol (float, optional) – Set tolerance for minimum distance between two nodes. The unit is the same as that of the watershed’s CRS. The default is 1.

  • meshpy.triangle.build() (Additional keyword arguments include all options for)

Condition#

watershed_workflow.condition.burnInRiver(river: River, network_burn_in_depth: Callable[[River], float]) → None[source]#

Reduce reach elevations by a float or function.

watershed_workflow.condition.computeChangeStatistics(original_vertex_elevations: ndarray, new_vertex_elevations: ndarray, metric: str = 'rmse') → float | Dict[str, float][source]#

Quantify deviation between original and modified vertex elevations.

Computes statistics describing how much vertex elevations changed, useful for assessing the impact of pit filling on the DEM.

Parameters:

  • original_vertex_elevations (np.ndarray) – Original vertex elevations, shape (num_vertices,).

  • new_vertex_elevations (np.ndarray) – Modified vertex elevations, shape (num_vertices,).

  • metric (str, optional) –

    Metric to compute. Options:

    • ’rmse’: Root mean square error (default)

    • ’mae’: Mean absolute error

    • ’max’: Maximum absolute change (worst case)

    • ’total’: Total (sum) of absolute changes

    • ’all’: Returns dict with all metrics plus additional statistics

    Default is ‘rmse’.

Returns:

If metric != ‘all’, returns single float value in mesh coordinate units.

If metric == ‘all’, returns dict with keys:

  • ’rmse’: Root mean square error, measures overall goodness of fit; penalizes large changes.

  • ’mae’: Mean absolute error across all vertices, a typical elevation change per vertex

  • ’max’: Maximum (worst-case) elevation change

  • ’total’: Sum of absolute elevation changes

  • ’num_modified’: Count of vertices with elevation change > 1e-9

  • ’mean_modified’: Mean change among modified vertices only

  • ’median_modified’: Median change among modified vertices only

  • ’percent_modified’: Percentage of vertices modified

Return type:

float or dict

Notes

All metrics are computed as (new - original), so positive values indicate vertices were raised. Pit filling should only raise elevations (or leave them unchanged), never lower them.

watershed_workflow.condition.conditionCell(m2: Mesh2D, c, pit: Tuple[int | bool, str, float, float], forced_outlet_edges: Set[Tuple[int, int]], optional_outlet_edges: Set[Tuple[int, int]], divide_edges: Set[Tuple[int, int]], epsilon: float, tol: float, additional_fixed_vertices: Set[int] | None = None, relative_to: Iterable[int] | None = None) → Dict[int, float][source]#

Compute vertex raises needed to condition a single pit cell.

Determines which vertices should be fixed based on boundary edge categorization, then computes appropriate raises to eliminate the pit.

Parameters:

  • m2 (Mesh2D) – The 2D mesh containing vertex coordinates and cell connectivity.

  • c (int) – Cell index into m2.conn.

  • pit (Tuple[int or bool, str, float, float]) – Pit information tuple: (is_pit, cause, internal_depth, boundary_depth).

  • forced_outlet_edges (Set[Tuple[int, int]]) – Boundary edges that must be outlets.

  • optional_outlet_edges (Set[Tuple[int, int]]) – Boundary edges that may be outlets.

  • divide_edges (Set[Tuple[int, int]]) – Boundary edges that must not be outlets (watershed divides).

  • epsilon (float) – Minimum elevation increase to enforce drainage.

  • tol (float) – Numerical tolerance for roundoff errors.

  • additional_fixed_vertices (Set[int] or None, optional) – Additional vertices that must remain fixed beyond those determined by boundary edge logic. Default is None (no additional fixed vertices).

  • relative_to (Iterable[int] or None, optional) – If provided, only consider these cells when measuring pits. Default is None (consider all neighbors).

Returns:

vertex_raises – Dictionary mapping vertex_id -> raise_amount. Returns empty dict if cell is already conditioned.

Return type:

Dict[int, float]

Notes

Logic by boundary edge type:

  • Internal cells: Raise centroid by internal_depth + epsilon

  • Forced outlet edges: Fix vertices on forced edges, raise centroid

  • Divide edges with boundary_depth > 0: Raise only boundary vertices

  • Divide edges with boundary_depth <= 0: Raise centroid by internal_depth

  • Optional outlet edges: Fix vertices on optional edges, raise centroid

Additional fixed vertices are respected in all cases.

watershed_workflow.condition.conditionMesh(m2: Mesh2D, preserved_pits: Iterable[int] | None = None, forced_outlet_edges: Iterable[Edge] | None = None, optional_outlet_edges: Iterable[Edge] | None = None, divide_edges: Iterable[Edge] | None = None, epsilon: float = 0.0, tol: float = 1e-08, plot: bool = False) → Tuple[Mesh2D, List[Dict[str, Any]]][source]#

The recommended algorithm for filling pits away from the river.

watershed_workflow.condition.conditionRiverMesh(m2: Mesh2D, river: River, smooth: bool = False, lower: bool = False, bank_integrity_elevation: float = 0.0, depress_headwaters_by: float | None = None, network_burn_in_depth: Callable[[River], float] | None = None, known_depressions: List[int] | None = None) → None[source]#

Condition, IN PLACE, the elevations of stream-corridor elements to ensure connectivity throgh culverts, skips ponds, maintain monotonicity, or otherwise enforce depths of constructed channels.

Parameters:

  • m2 (watershed_workflow.mesh.Mesh2D object) – 2D mesh with 3D coordinates.

  • river (watershed_workflow.river_tree.River object) – River tree with reach[‘elems’] added for quads

  • smooth (boolean, optional) – If true, smooth the profile of each reach using a gaussian filter (mainly to pass through railroads and avoid reservoirs).

  • lower (boolean, optional) – If true, lower the smoothed bed profile to match the lower points on the raw bed profile. This is useful particularly for narrow ag. ditches where NHDPLus flowlines often do not coincide with the DEM depressions and so stream-elements intermitently fall into them.

  • bank_integrity_elevation (float, optional) – Where the river is passing right next to the reservoir or NHDline is misplaced into the reservoir, banks may fall into the reservoir. If true, this will enforce that the bank vertex is at a higher elevation than the stream bed elevation.

  • depress_headwaters_by (float, optional) – If the depression is not captured well in the DEM, the river-mesh elements (streambed) headwater reaches may be lowered by this number. The effect is propogated downstream only up to where it is needed to maintain topographic gradients on the network scale in the network sweep step.

  • network_burn_in_depth (Callable[[River,], float], optional) – A function that takes a reach (River object) as input and returns the burn-in depth for that specific reach. This depth specifies how much to lower the river-mesh elements below their original elevation. The callable allows for dynamic calculation based on reach properties, stream order, or custom logic.

  • known_depressions (list, optional) – If provided, a list of IDs to not be burned in via the network sweep.

watershed_workflow.condition.conditionRiverMeshes(m2: Mesh2D, rivers: List[River], *args, **kwargs) → None[source]#

For multiple rivers, condition, IN PLACE, the elevations of stream-corridor elements to ensure connectivity throgh culverts, skips ponds, maintain monotonicity, or otherwise enforce depths of constructed channels.

watershed_workflow.condition.distributeProfileToMesh(m2: Mesh2D, river: River) → None[source]#

Take reach profile elevations and move them out to the mesh vertices.

watershed_workflow.condition.enforceBankIntegrity(m2: Mesh2D, river: River, bank_integrity_elevation: float) → None[source]#

Forces banks at least bank_integrity_elevation higher than the channel elevation.

watershed_workflow.condition.enforceLocalMonotonicity(reach: River, moving: Literal['downstream', 'upstream'] = 'downstream') → None[source]#

Ensures that the streambed-profile elevations are monotonically increasing as we move upstream, or decreasing as we move downstream.

watershed_workflow.condition.enforceMonotonicity(river: River, depress_headwaters_by: float | None = None, known_depressions: List[int] | None = None) → None[source]#

Sweep the river network from each headwater reach (leaf node) to the watershed outlet (root node), removing aritificial obstructions in the river mesh and enforcing depths of constructed channels.

watershed_workflow.condition.fillPits(m2: Mesh2D, method_name: str = 'marching', preserved_pits: Iterable[int] | None = None, forced_outlet_edges: Iterable[Edge] | None = None, optional_outlet_edges: Iterable[Edge] | None = None, divide_edges: Iterable[Edge] | None = None, epsilon: float = 0.0, tol: float = 1e-08, plot: bool = False, max_iterations: int | None = None, **kwargs) → Tuple[Mesh2D, Dict[str, Any]][source]#

Fill pits in mesh using specified method(s) with optional plotting and metrics.

User-friendly wrapper that runs pit filling method(s), computes metrics, and optionally plots results. Note that the returned Mesh2D may be the input mesh, or may be different, depending upon the algorithm. The user should assume that m2 is modified in place, but should use the returned mesh.

Parameters:

  • m2 (Mesh2D) – The mesh to condition (modified in place).

  • method_name (str, optional) – Method to use. Can be: ‘recommended’, ‘global’, ‘marching’, ‘marching old’, ‘null’, or ‘boundary cleanup’. Default is ‘recommended’.

  • preserved_pits (Iterable[int], optional) – Cell indices to preserve as pits (e.g., lakes). Default is the empty list.

  • forced_outlet_edges (Iterable[Tuple[int, int]], optional) – Boundary edges that must be outlets. Default is the empty list.

  • optional_outlet_edges (Iterable[Tuple[int, int]], optional) – Boundary edges that may be outlets. Defaults to all boundary edges not in forced_outlet_edges or divide_edges.

  • divide_edges (Iterable[Tuple[int, int]], optional) – Boundary edges that must not be outlets (watershed divides). Default is the empty list.

  • epsilon (float, optional) – Minimum slope parameter. Default is 0.0 (no enforced slope).

  • tol (float, optional) – Numerical tolerance for roundoff errors. Default is 1.e-8.

  • plot (bool, optional) – If True, creates before/after elevation comparison plots. Default is False.

  • max_iterations (int, optional) – Maximum iterations for iterative methods. If specified, uses iterative version of the method. Default is None (use non-iterative version).

  • kwargs (dict, optional) – Additional keyword arguments passed to fillPits algorithm.

Returns:

  • m2 (Mesh2D) – The conditioned mesh.

  • result (Dict[str, Any]) – Statistics dictionary with keys: - ‘method_name’: Name of method used - ‘pits_initial’: List of initial pits - ‘pits_final’: List of remaining pits - ‘pits_removed’: Number of pits removed - ‘elevation_stats’: Dict with ‘rmse’, ‘mae’, ‘max’, ‘num_modified’, etc.

watershed_workflow.condition.fillPits_global(m2: Mesh2D, pits: List[Tuple[int, str, float, float]], preserved_pits: Set[int], forced_outlet_edges: Set[Edge], optional_outlet_edges: Set[Edge], divide_edges: Set[Edge], epsilon: float, tol: float, boundary_only: bool = False) → Mesh2D[source]#

Fill pits in the mesh using iterative cell-based algorithm.

This algorithm iteratively identifies and eliminates pits by raising vertex elevations using boundary edge categorization to determine which vertices should be fixed. Guarantees no pits remain (within tolerance) and only raises elevations, never lowering them.

Modifies m2.coords[:, 2] in place.

Parameters:

  • m2 (Mesh2D) – The 2D mesh containing vertex coordinates and cell connectivity. Elevations in m2.coords[:, 2] will be modified in place.

  • pits (List[Tuple[int, str, float, float]]) – List of pits to fill: (cell, cause, internal_depth, boundary_depth).

  • preserved_pits (Set[int]) – Cell indices that should be preserved as pits (e.g., lakes, playas, or other real depressions). These cells are never filled.

  • forced_outlet_edges (Set[Tuple[int, int]]) – Boundary edges that must be outlets. Cells touching these edges will have these edge vertices fixed while raising the cell centroid.

  • optional_outlet_edges (Set[Tuple[int, int]]) – Boundary edges that may be outlets.

  • divide_edges (Set[Tuple[int, int]]) – Boundary edges that must not be outlets (watershed divides). Cells touching these may have divide edge vertices raised to prevent outward flow.

  • epsilon (float) – Minimum elevation increase per cell in flow direction. Units are same as mesh coordinates.

  • tol (float) – Numerical tolerance for roundoff errors.

  • boundary_only (bool, optional) – If True, only process boundary pits (pits with cause != ‘internal’). If False, process all pits. Default is False.

Notes

Algorithm iteratively:

  1. Identifies pit cells using findPits() with boundary edge categorization

  2. For each pit, calls computePitDepth() to get internal and boundary depths

  3. Calls conditionCell() to determine which vertices to raise based on boundary edge type (forced outlet, divide, or optional outlet)

  4. Accumulates vertex raises (taking max for vertices shared by multiple cells)

  5. Applies vertex elevation changes to mesh

  6. Repeats until no changes exceed tolerance or max_iterations reached

Boundary edge categorization ensures: - Forced outlet edges: Interior vertices raised, boundary vertices fixed - Divide edges: Boundary vertices raised to prevent outward flow - Optional outlet edges: Interior vertices raised, boundary vertices fixed

Elevations are only raised, never lowered.

watershed_workflow.condition.fillPits_marching(m2: Mesh2D, pits: List[Tuple[int, str, float, float]], preserved_pits: Set[int], forced_outlet_edges: Set[Edge], optional_outlet_edges: Set[Edge], divide_edges: Set[Edge], epsilon: float, tol: float, seed_policy: str | List[str] | None = None, seed_to: str | None = 'waterway', conditioning_policy: str = 'waterway', fixing_policy: str = 'waterway', preserved_pits_are_fixed: bool = True, replace_upon_conditioning: bool = True) → Mesh2D[source]#

Fill pits using a greedy marching algorithm.

The goal of this algorithm is to ensure that, starting from an outlet cell and a list of known pits, there is a path to every cell by way of faces that are monotonically increasing in elevation.

A cell is called reachable if such a path exists. Cells are incrementally added to the “waterway,” or the set of cells with identified, unchanging paths.

This algorithm starts with all preserved pits and outlets. These cells are fixed and placed in the waterway. Cells are added to the waterway by picking the lowest elevation cell that currently borders the existing waterway.

A cell must be conditioned before it can be added to the waterway. Conditioning a cell enforces that this cell has a higher elevation than at least one of its neighbors that is already in the waterway.

If a cell has a lower cell that is NOT in the waterway, and has a valid, monotonically decreasing pathway to the waterway, then that pathway has a cell that is in the border of the current waterway, and that cell’s elevation is lower than this cell, and therefore would have been selected before this cell. Contradiction; therefore the downhill path must lead to a pit or the boundary, but not to the waterway. This observation is key to enforcing the condition.

To ensure that conditioning one cell does not break other, already conditioned cells, we must choose a set of vertices to fix. Clearly all cells in the waterway should have their vertices fixed; these are not revisited. The obvious approach is to fix vertices upon conditioning – then once a cell is conditioned, it need not be revisited.

The only known failure mechanism of this approach is that a cell cannot be conditioned because the vertices of that cell are all a part of previously conditioned cells and therefore have been fixed. In that case the cell is added anyway and left unconditioned (becoming a pit).

Note there are lot of options for the algorithm, but the defaults for all policies are expected to be the most robust.

Parameters:

  • m2 (Mesh2D) – The mesh to condition.

  • pits (List[Tuple[int, str, float, float]]) – List of pits to fill: (cell, cause, internal_depth, boundary_depth).

  • preserved_pits (Set[int]) – Cell indices to preserve as pits (e.g., lakes) that may be a depression. This is important for reservoirs/lakes/etc where bathymetry is known and pits are physical.

  • forced_outlet_edges (Set[Tuple[int, int]]) – Boundary edges that are forced outlets.

  • optional_outlet_edges (Set[Tuple[int, int]]) – Boundary edges that may be outlets.

  • divide_edges (Set[Tuple[int, int]]) – Boundary edges that are watershed divides.

  • epsilon (float) – Minimum slope parameter.

  • tol (float) – Tolerance for numeric roundoff.

  • seed_policy (str or List[str], optional) – What elements are included in the initial seed? Default is [‘forced outlets’, ‘preserved pits’, ‘optional outlets’]. Valid entries include these values.

  • seed_to (str, optional) – Where are seeded elements put, into the ‘waterway’ or into the ‘border’. Default is ‘waterway’.

  • conditioning_policy (str, optional) – When to condition – upon entering ‘waterway’ or ‘border’. Default is ‘waterway’.

  • fixing_policy (str, optional) – When to fix vertices – upon entering ‘waterway’ or ‘border’. Default is ‘waterway’.

  • preserved_pits_are_fixed (bool, optional) – If True, fixes all vertices of preserved pits. Default is True.

  • replace_upon_conditioning (bool, optional) – If cells are conditioned on placement in the waterway, conditioning may raise the cell elevation, meaning it is no longer the lowest elevation cell in the border. If this is True, conditioned cells are placed back into the border and not put in the waterway. Default is True.

watershed_workflow.condition.fillPits_marching_old(m2: Mesh2D, pits: List[Tuple[int, str, float, float]], preserved_pits: Set[int], forced_outlet_edges: Set[Edge], optional_outlet_edges: Set[Edge], divide_edges: Set[Edge], epsilon: float, tol: float) → Mesh2D[source]#

Fill pits using a greedy marching algorithm.

The goal of this algorithm is to ensure that, starting with an outlet cell and a list of known pits, there is a path to every cell by way of faces that is monotonically increasing in elevation.

A cell is called reachable if such a path exists. Cells are incrementally added to the “waterway,” or the set of cells with identified paths.

Starting from an outlet, it adds cells to the waterway by picking the lowest elevation cell that currently borders the existing waterway. It conditions upon adding the cell to the boundary.

Conditioning a cell requires that this cell is higher that at least one of its neighbors that is already in the waterway. Note that this is more aggressive because it requires the lower cell to be in the waterway, not just any neighbor.

If a cell has a lower cell that is NOT in the waterway, and has a valid, monotonically decreasing pathway to the waterway, then that pathway has a cell that is in the boundary of the current waterway, and that cell’s elevation is lower than this cell, and therefore would have been selected before this cell. Contradiction; therefore the downhill path must lead to a pit or the boundary, but not to the waterway.

Parameters:

  • m2 (Mesh2D) – The mesh to condition.

  • pits (List[Tuple[int, str, float, float]]) – List of pits to fill: (cell, cause, internal_depth, boundary_depth).

  • preserved_pits (Set[int]) – Cell indices to preserve as pits (e.g., lakes) that may be a depression. This is important for reservoirs/lakes/etc where bathymetry is known and pits are physical.

  • forced_outlet_edges (Set[Tuple[int, int]]) – Boundary edges that are forced outlets.

  • optional_outlet_edges (Set[Tuple[int, int]]) – Boundary edges that may be outlets.

  • divide_edges (Set[Tuple[int, int]]) – Boundary edges that are watershed divides.

  • epsilon (float) – Minimum slope parameter.

  • tol (float) – Numerical tolerance for roundoff errors.

watershed_workflow.condition.findPits(m2: Mesh2D, preserved_pits: Iterable[int] | None = None, forced_outlet_edges: Iterable[Tuple[int, int]] | None = None, optional_outlet_edges: Iterable[Tuple[int, int]] | None = None, divide_edges: Iterable[Tuple[int, int]] | None = None, epsilon: float = 0.0, tol: float = 1e-08) → List[Tuple[int, str, float, float]][source]#

Identify problematic pits (local minima) in the mesh.

Finds cells whose centroid elevation is lower than their surroundings, preventing drainage. Uses sophisticated boundary edge categorization to handle outlet and divide edges correctly.

Parameters:

  • m2 (Mesh2D) – The 2D mesh containing vertex coordinates and cell connectivity. Elevations are read from m2.coords[:, 2].

  • preserved_pits (Iterable[int], optional) – Cell indices that are intentionally pits (e.g., lakes) and should not be reported as problems. Default is the empty list.

  • forced_outlet_edges (Iterable[Tuple[int, int]], optional) – Boundary edges that must be outlets. Cells touching these edges are pits only if water would flow inward across the boundary. Default is the empty list.

  • optional_outlet_edges (Iterable[Tuple[int, int]], optional) – Boundary edges that may be outlets. Defaults to all boundary edges not in forced_outlet_edges or divide_edges. Cells touching these are pits only if trapped both internally and externally.

  • divide_edges (Iterable[Tuple[int, int]], optional) – Boundary edges that must not be outlets (watershed divides). Cells touching these are pits if water would flow outward OR if trapped internally. Default is the empty list.

  • epsilon (float, optional) – Minimum elevation increase required to not be a pit. Use the same epsilon as fillPits methods for consistency. Default is 0.0.

  • tol (float, optional) – Numerical tolerance for roundoff errors. Default is 1.e-8.

Returns:

pits – List of pits: [cell, cause, internal_pit_depth, boundary_pit_depth]

Return type:

List[Tuple[int, str, float, float]]

Notes

Pit detection uses: isPit(depth) = depth > -(epsilon - tol) where depth is from computePitDepth().

Boundary edge partitioning:

  • All boundary edges are categorized into forced_outlet_edges, optional_outlet_edges, or divide_edges

  • Unspecified edges default to optional_outlet_edges

  • forced_outlet_edges take precedence in overlaps

Pit criteria by cell type:

  • Internal cells: isPit(internal_depth)

  • Forced outlet cells: isPit(boundary_depth)

  • Divide edge cells: isPit(internal_depth) OR NOT isPit(boundary_depth)

  • Optional outlet cells: isPit(internal_depth) AND isPit(boundary_depth)

watershed_workflow.condition.plotPitFilling(m2: Mesh2D, old_pits: List[Tuple[int, str, float, float]], new_pits: List[Tuple[int, str, float, float]], old_verts: ndarray, method_name: str, ax: Any | None = None, metrics: Dict | None = None, pit_kwargs=None, vertex_kwargs=None, cell_kwargs=None)[source]#

Plot before/after elevation comparison for pit filling algorithms.

Creates 3x3 panel plot showing: - Row 0: Pit depths (old, new, delta) - Row 1: Vertex elevations (old, new, delta) - Row 2: Cell centroid elevations (old, new, delta)

Includes dynamic colormap scaling that automatically adjusts color limits when zooming/panning the plot for detailed inspection of specific regions.

Parameters:

  • m2 (Mesh2D) – The mesh with new elevations in m2.coords[:,2]

  • old_pits (List[Tuple[int, str, float, float]]) – List of pits before filling

  • new_pits (List[Tuple[int, str, float, float]]) – List of pits after filling

  • old_verts (np.ndarray) – Copy of original vertex before modification

  • method_name (str) – Name of the algorithm (for title)

  • ax (array of matplotlib axes, optional) – If provided, should be 3x3 array of axes

  • metrics (dict, optional) – Dictionary with metrics to display on plot. Expected keys: ‘pits_initial’, ‘pits_final’, ‘rmse’, ‘mae’, ‘max’, ‘num_modified’

  • pit_kwargs (dict, optional) – Additional arguments passed to m2.plot() for pit depth plots

  • vertex_kwargs (dict, optional) – Additional arguments passed to m2.plotVertices() for vertex plots

  • cell_kwargs (dict, optional) – Additional arguments passed to m2.plot() for cell plots

Returns:

  • fig (matplotlib.figure.Figure) – The figure object

  • axes (np.ndarray) – Array of axes (3x3)

Notes

Dynamic colormap scaling groups related plots: - Pit depths (old/new) share linear scale, delta uses symmetric scale - Elevations (old/new vertex and cell) share linear scale - Delta elevations use symmetric scales around zero

watershed_workflow.condition.raiseCellCentroid(m2: Mesh2D, c: int, target_raise: float, fixed_vertices: Set[int]) → Dict[int, float][source]#

Compute vertex raises to increase a cell’s centroid by target amount.

Raises only non-fixed vertices uniformly so that the cell’s centroid elevation increases by exactly target_raise. The centroid is computed as the arithmetic mean of vertex coordinates.

Parameters:

  • m2 (Mesh2D) – The 2D mesh containing vertex coordinates and cell connectivity.

  • c (int) – Cell index into m2.conn.

  • target_raise (float) – Amount to raise the cell centroid elevation.

  • fixed_vertices (Set[int]) – Set of vertex indices that cannot be raised (typically vertices on boundary edges that must remain at their current elevation).

Returns:

vertex_raises – Dictionary mapping vertex_id -> raise_amount for free vertices only. Returns empty dict if no free vertices exist (all vertices fixed).

Return type:

Dict[int, float]

Notes

To raise centroid by target_raise with n_fixed fixed vertices: raise = target_raise * n_total / n_free

This ensures that the mean of all vertex elevations increases by exactly target_raise when only the free vertices are modified.

watershed_workflow.condition.setProfileByDEM(rivers: List[River], dem: DataArray, **kwargs) → None[source]#

Set the z-coordinate of the reach linestring from a DEM dataset.

watershed_workflow.condition.singlePitDepth(p)[source]#

Given a pit-tuple, returns a single depth used for debugging and info.

watershed_workflow.condition.smoothProfile(reach: River, lower: bool = False) → None[source]#

Applies gaussian filter smoothing to the bed-profile obtained from DEM.

This option becomes important in ag. watersheds when NHDPLus is inconsistent with the depression in the DEM.

Regions#

This namespace hosts functions to add labelsets to define regions

watershed_workflow.regions.addDischargeRegions(m2, discharge_points, labels=None, include_cells=True, buffer_width=1)[source]#

Add labeled sets for three faces for each discharge point in the river corridor. The three faces include downstream shorter edge of the quad and two edges connecting two downstream vertices of the quad and non-quad vertice on the bank triangle. Corresponding upstreams cells (a quad and two triangles) are also added if include_cells is True, which should be use with <Parameter name=”direction normalized flux relative to region” type=”string” value=”discharge_cell_region_name” /> in the ATS observation parameter list in the input file

Parameters:

  • m2 (watershed_workflow.mesh.Mesh2D) – The 2D mesh containing river corridor elements

  • discharge_points (list of (x,y) coordinates) – List of discharge point locations to add regions for

  • labels (list of str, optional) – Custom labels for each discharge point. If not provided, defaults to ‘discharge point 0’, ‘discharge point 1’, etc.

  • include_cells (bool, optional) – If True, add a labeled set for the cells just upstream of discharge faces. Default is True.

  • buffer_width (float, optional) – Buffer width to identify quad elements containing the discharge point. Default is 1.

Notes

For each discharge point, this creates a labeled set containing three edges: 1. The downstream edge of the quad element containing the point 2. Edge connecting downstream right vertex to right bank 3. Edge connecting downstream left vertex to left bank and labeled set containing three cells: 1. the quad element containing the point 2. the two triangles sharing edges with the quad

watershed_workflow.regions.addPolygonalRegions(m2: Mesh2D, polygons: SplitHUCs | List[Tuple[Polygon, str]], volume: bool = False, return_partitions: bool = False) → List[List[int]] | None[source]#

Label m2 with region(s) for each polygon.

Always adds a surface region; if volume also adds a volume region for extrusion in 3D.

Parameters:

  • m2 (mesh.Mesh2D) – The mesh to label.

  • polygons (SplitHUCs | List[Tuple[shapely.Polygon, str]]) – The polygons covering each region. Either a SplitHUCs object or a list of tuples containing (polygon, name) pairs.

  • volume (bool, optional) – If true, also add the volumetric region below the polygon that will be extruded in a 3D mesh eventually.

  • return_partitions (bool, optional) – If true, return the list of cell indices for each polygon. Default is False.

Returns:

partitions – A list of length polygons, each entry of which is the list of cell indices in that polygon. Only returned if return_partitions is True.

Return type:

List[List[int]] or None

watershed_workflow.regions.addReachIDRegions(m2: Mesh2D, river: River)[source]#

Add labeled sets to m2 for reaches of each stream order .

Parameters:#

m2: Mesh2D

2D mesh elevated on DEMs

river: watershed_workflow.river_tree.RiverTree

List of rivers used to create the river corridors.

reaches: list(str)

list of NHDID IDs to be labeled.

watershed_workflow.regions.addRiverCorridorRegions(m2: Mesh2D, rivers: List[River], labels: List[str] | None = None)[source]#

Add labeled sets to m2 for each river corridor.

Parameters:#

m2: Mesh2D

2D mesh elevated on DEMs

rivers: list(watershed_workflow.river_tree.RiverTree)

List of rivers used to create the river corridors.

labels: list(str), optional

List of names, one per river.

watershed_workflow.regions.addStreamOrderRegions(m2: Mesh2D, rivers: List[River]) → None[source]#

Add labeled sets to m2 for reaches of each stream order .

Parameters:#

m2: Mesh2D

2D mesh elevated on DEMs

riversList[watershed_workflow.river_tree.River]

River used to create the river corridors.

watershed_workflow.regions.addSurfaceRegions(m2: Mesh2D, column: str = 'land_cover', names: dict[int, str] | None = None)[source]#

Add labeled sets to a mesh – one per unique color.

Parameters:

  • m2 (mesh.Mesh2D) – The mesh to label.

  • names (dict, optional) – Dictionary mapping colors to color name.

watershed_workflow.regions.addWatershedAndOutletRegions(m2: Mesh2D, hucs: SplitHUCs, outlet_width: float = 300, exterior_outlet=True)[source]#

Add four labeled sets to m2 for each polygon:

  • cells in the polygon, to be extruded

  • cells in the polygon, to be kept as faces upon extrusion

  • boundary of the polygon (edges, kept as faces)

  • outlet of the polygon (edges, kept as faces, within outlet width of the outlet)

Parameters:

  • m2 (mesh.Mesh2D) – The mesh to label.

  • hucs (iterable[shapely.Polygon] or SplitHUCs) – Watershed polygons.

  • outlet_width (float, optional) – How wide should the outlet region be? Note this should include not just the river width, but more realistically the best resolved floodplain, or 1-2 face-widths, whichever is bigger.

  • exterior_outlet (bool, optional) – If true, find the outlet point that intersects the boundary and include regions around that outlet as well.

watershed_workflow.regions.findDischargeEdgesCells(m2, discharge_point, include_cells=True, buffer_width=1)[source]#

Find the edges around a discharge point in a river corridor mesh.

Parameters:

  • m2 (watershed_workflow.mesh.Mesh2D) – The 2D mesh containing river corridor elements

  • discharge_point (shapely.geometry.Point) – The discharge point location

  • include_cells (bool, optional) – If True, include cells in the labeled set. Default is True.

  • buffer_width (float, optional) – Buffer width to identify quad elements containing the discharge point. Default is 1.

Returns:

List of (vertex1, vertex2) pairs defining edges. if include_cells is True, also returns the cells just upstream of the edges.

Return type:

list of tuples