Source code for watershed_workflow.triangulation

"""Triangulates polygons"""
from typing import Tuple, List, Optional, Set, Iterator

import logging
import collections
import itertools
import numpy as np
import numpy.linalg as la
from matplotlib import pyplot as plt
import scipy.spatial

import shapely

import watershed_workflow.river_tree
import watershed_workflow.split_hucs


[docs] class Nodes: """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. """ def __init__(self, decimals: int = 3) -> None: """Initialize a Nodes collection. Parameters ---------- decimals : int, optional Number of decimal places for rounding coordinates. Default is 3. """ self.decimals = decimals self._i = 0 self._store: collections.OrderedDict[Tuple[float, ...], int] = collections.OrderedDict() def __len__(self) -> int: """Length is tracked by the counter""" return self._i def __getitem__(self, key: Tuple[float, ...]) -> int: """Get item index, setting if default.""" # Note this implementation could be optimized, and should be # done so assuming keys are NOT there, which is the most # likely case. key = tuple(round(p, self.decimals) for p in key) if key not in self._store: self._store[key] = self._i self._i += 1 return self._store[key] def __iter__(self) -> Iterator[Tuple[float, ...]]: """Iterable collection""" for k in self._store.keys(): yield k
[docs] def connectOnewayTrip(inds: List[int]) -> List[Tuple[int, int]]: """Connect indices in edges in a oneway fashion""" return [(inds[i], inds[i + 1]) for i in range(len(inds) - 1)]
[docs] def connectRoundTrip(inds: List[int]) -> List[Tuple[int, int]]: """Connect indices in edges in a round-trip fashion""" return connectOnewayTrip(inds) + [(inds[-1], inds[0]), ]
[docs] def orient(e: Tuple[int, int]) -> Optional[Tuple[int, int]]: """Orient an edge consistently. Parameters ---------- e : Tuple[int, int] Edge as tuple of two vertex indices. Returns ------- Optional[Tuple[int, int]] Oriented edge or None if vertices are identical. """ if e[0] > e[1]: return e[1], e[0] elif e[0] < e[1]: return e[0], e[1] else: return None
[docs] class NodesEdges: """A collection of nodes and edges.""" def __init__( self, objlist: Optional[List[shapely.geometry.LineString | shapely.geometry.Polygon]] = None ) -> None: """Initialize NodesEdges collection. Parameters ---------- objlist : Optional[List[shapely.geometry.LineString | shapely.geometry.Polygon]], optional List of geometry objects to add. Default is None. """ self.nodes = Nodes() self.edges: Set[Tuple[int, int]] = set() if objlist != None: for obj in objlist: self.add(obj)
[docs] def add(self, obj: shapely.geometry.LineString | shapely.geometry.Polygon) -> None: """Adds nodes and edges from obj into collection.""" if isinstance(obj, shapely.geometry.LineString): inds = [self.nodes[c] for c in obj.coords] for e in connectOnewayTrip(inds): oe = orient(e) if oe is not None: self.edges.add(oe) elif isinstance(obj, shapely.geometry.Polygon): inds = [self.nodes[c] for c in obj.boundary.coords] for e in connectRoundTrip(inds): oe = orient(e) if oe is not None: self.edges.add(oe) else: raise TypeError("Invalid type for add, %r" % type(obj))
[docs] def check(self, tol: float) -> None: """Checks consistency of the internal representation.""" logging.info(" checking graph consistency") logging.info(" tolerance is set to {}".format(tol)) min_dist = 1.e10 coords = np.array(list(self.nodes)) kdtree = scipy.spatial.cKDTree(coords) bad_pairs = kdtree.query_pairs(tol) # Retrieve the coordinates for each bad pair bad_pair_coords = [(coords[i], coords[j]) for i, j in bad_pairs] if len(bad_pairs) != 0: raise ValueError('tol= {} is too large, try decrease tolerance!'.format(tol) + 'or check bad pairs={}'.format(bad_pair_coords)) min_node = min(self.nodes[n] for n in self.nodes) max_node = max(self.nodes[n] for n in self.nodes) assert (min_node == 0) assert (max_node == len(self.nodes) - 1) min_edge_node = min(n for e in self.edges for n in e) max_edge_node = max(n for e in self.edges for n in e) assert (min_edge_node == 0) assert (max_edge_node == len(self.nodes) - 1)
[docs] def triangulate(hucs: watershed_workflow.split_hucs.SplitHUCs, internal_boundaries: Optional[List[shapely.geometry.LineString]] = None, hole_points: Optional[List[shapely.geometry.Point]] = None, tol: float = 1.0, **kwargs) -> Tuple[np.ndarray, np.ndarray]: """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). 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. Additional keyword arguments include all options for meshpy.triangle.build() """ import meshpy.triangle logging.info("Triangulating...") # get a list of required boundaries, both interior and exterior linestrings: List[shapely.geometry.LineStrings] = [] # -- internal boundaries may include a river mesh and must go first if internal_boundaries is not None: linestrings.extend(internal_boundaries) linestrings = linestrings + list(hucs.linestrings) # convert to nodes and edges nodes_edges = NodesEdges(linestrings) logging.info(" %i points and %i facets" % (len(nodes_edges.nodes), len(nodes_edges.edges))) nodes_edges.check(tol=tol) logging.info(" building graph data structures") info = meshpy.triangle.MeshInfo() nodes = np.array(list(nodes_edges.nodes)) pdata = [tuple([float(c) for c in p]) for p in nodes] info.set_points(pdata) fdata = [[int(i) for i in f] for f in nodes_edges.edges] info.set_facets(fdata) # add hole points, which should include the river mesh interior if hole_points is not None: info.set_holes(hole_points) logging.info(" triangle.build...") # pop this option if false, which silences the warning if it does # not exist but we didn't ask for it anyway. if 'enforce_delaunay' in kwargs.keys() and not kwargs['enforce_delaunay']: kwargs.pop('enforce_delaunay') try: mesh = meshpy.triangle.build(info, **kwargs) except TypeError as err: try: # our modification to meshpy.triangle is not present, try without it kwargs.pop('enforce_delaunay') except KeyError: raise err else: logging.warning( "Triangulate: '--enforce-delaunay' option requires a hacked `meshpy.triangle`. Proceeding without this option because it is not recognized." ) mesh = meshpy.triangle.build(info, **kwargs) mesh_points = np.array(mesh.points) mesh_tris = np.array(mesh.elements) logging.info(" ...built: %i mesh points and %i triangles" % (len(mesh_points), len(mesh_tris))) return mesh_points, mesh_tris
[docs] def refineByMaxArea(max_area): """Returns a refinement function based on max area, for use with Triangle.""" def refine(vertices, area): """A function for use with watershed_workflow.triangulate.triangulate's refinement_func argument based on a global max area.""" res = bool(area > max_area) # if area < 1.e-5: # raise RuntimeError("TinyTriangle Error") return res return refine
[docs] def refineByRiverDistance(near_distance, near_area, away_distance, away_area, rivers): """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. """ def max_area_valid(distance): """A function to make sure max area scales with distance from river network Units in [m] """ if distance > away_distance: area = away_area elif distance < near_distance: area = near_area else: area = near_area + (distance-near_distance) / (away_distance -near_distance) * (away_area-near_area) return area if type(rivers[0]) == shapely.geometry.Polygon: river_multiline = shapely.geometry.MultiPolygon(rivers) else: river_multiline = shapely.geometry.MultiLineString( [r.linestring for river in rivers for r in river]) def refine(vertices, area): """A function for use with watershed_workflow.triangulate.triangulate's refinement_func argument based on size gradation from a river.""" bary = np.sum(np.array(vertices), axis=0) / 3 bary_p = shapely.geometry.Point(bary[0], bary[1]) distance = bary_p.distance(river_multiline) max_area = max_area_valid(distance) res = bool(area > max_area_valid(distance)) #logging.debug("refine? area = %g, distance = %g, max_area = %g: refine = %r"%(area,distance,max_area,res)) return res return refine
[docs] def refineByPolygons(polygons, areas): """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. """ def refine(vertices, area): """A function for use with watershed_workflow.triangulate.triangulate's refinement_func argument based on polygon area limits.""" bary = np.sum(np.array(vertices), axis=0) / 3 bary_p = shapely.geometry.Point(bary[0], bary[1]) # Check if a single area value is provided if isinstance(areas, (int, float)): max_area = areas for polygon in polygons: if polygon.contains(bary_p) and area > max_area: return True else: # Assume areas is a list of area limits for polygon, max_area in zip(polygons, areas): if polygon.contains(bary_p) and area > max_area: return True return False return refine
[docs] def refineByMaxEdgeLength(edge_length): """Returns a refinement function based on max edge length, for use with Triangle.""" def refine(vertices, area): verts4 = np.array([vertices[0], vertices[1], vertices[2], vertices[0]]) edge_lengths = la.norm(verts4[1:] - verts4[:-1], 2, 1) return bool(edge_lengths.max() > edge_length) return refine