Geometry

DHyDAMO Geometry API

common

as_linestring_list(linestring: LineString | MultiLineString | list[LineString | MultiLineString]) -> list[LineString]

Returns a list of LineStrings from a given LineString or MultiLineString. Useful for iterating over varying multi of single type geometries.

Parameters:

Name	Type	Description	Default
linestring	Union[LineString, MultiLineString]	LineString of MultiLineString shapely geometry object	required

Returns:

Type	Description
list[LineString]	List[LineString]: list of shapely geometry LineString objects

Source code in hydrolib/dhydamo/geometry/common.py

def as_linestring_list(
    linestring: LineString | MultiLineString | list[LineString | MultiLineString]
) -> list[LineString]:
    """Returns a list of LineStrings from a given LineString or MultiLineString. Useful for
    iterating over varying multi of single type geometries.

    Args:
        linestring (Union[LineString, MultiLineString]): LineString of MultiLineString shapely geometry object

    Returns:
        List[LineString]: list of shapely geometry LineString objects
    """
    return _as_geometry_list(linestring, LineString, MultiLineString)

as_point_list(point: Point | MultiPoint | list[Point | MultiPoint]) -> list[Point]

Returns a list of Point from a given Point or MultiPoint. Useful for iterating over varying multi of single type geometries.

Parameters:

Name	Type	Description	Default
linestring	Union[Point, MultiPoint]	Point of MultiPoint shapely geometry object	required

Returns:

Type	Description
list[Point]	List[Point]: list of shapely geometry Point objects

Source code in hydrolib/dhydamo/geometry/common.py

def as_point_list(
    point: Point | MultiPoint | list[Point | MultiPoint]
) -> list[Point]:
    """Returns a list of Point from a given Point or MultiPoint. Useful for
    iterating over varying multi of single type geometries.

    Args:
        linestring (Union[Point, MultiPoint]): Point of MultiPoint shapely geometry object

    Returns:
        List[Point]: list of shapely geometry Point objects
    """
    return _as_geometry_list(point, Point, MultiPoint)

as_polygon_list(polygon: Polygon | MultiPolygon | list[Polygon | MultiPolygon]) -> list[Polygon]

Returns a list of Polygons from a given Polygon or MultiPolygon. Useful for iterating over varying multi of single type geometries.

Parameters:

Name	Type	Description	Default
linestring	Union[Polygon, MultiPolygon]	Polygon of MultiPolygon shapely geometry object	required

Returns:

Type	Description
list[Polygon]	List[Polygon]: list of shapely geometry Polygon objects

Source code in hydrolib/dhydamo/geometry/common.py

def as_polygon_list(
    polygon: Polygon | MultiPolygon | list[Polygon | MultiPolygon]
) -> list[Polygon]:
    """Returns a list of Polygons from a given Polygon or MultiPolygon. Useful for
    iterating over varying multi of single type geometries.

    Args:
        linestring (Union[Polygon, MultiPolygon]): Polygon of MultiPolygon shapely geometry object

    Returns:
        List[Polygon]: list of shapely geometry Polygon objects
    """
    return _as_geometry_list(polygon, Polygon, MultiPolygon)

mesh

links1d2d_add_links_1d_to_2d(network: Network, branchids: list[str] = None, within: Polygon | MultiPolygon = None, max_length: float = np.inf) -> None

Function to add 1d2d links to network, by generating them from 1d to 2d. Branchids can be specified for 1d branches that need to be linked. A (Multi)Polygon can be provided were links should be made.

Parameters:

Name	Type	Description	Default
network	Network	Network in which the connections are made	required
branchids	List[str]	List of branchid's to connect. If None, all branches are connected. Defaults to None.	None
within	Union[Polygon, MultiPolygon]	Area within which connections are made. Defaults to None.	None
max_length	float	Max edge length. Defaults to None.	inf

Source code in hydrolib/dhydamo/geometry/mesh.py

def links1d2d_add_links_1d_to_2d(
    network: Network,
    branchids: list[str] = None,
    within: Polygon | MultiPolygon = None,
    max_length: float = np.inf,
) -> None:
    """Function to add 1d2d links to network, by generating them from 1d to 2d.
    Branchids can be specified for 1d branches that need to be linked.
    A (Multi)Polygon can be provided were links should be made.

    Args:
        network (Network): Network in which the connections are made
        branchids (List[str], optional): List of branchid's to connect. If None, all branches are connected. Defaults to None.
        within (Union[Polygon, MultiPolygon], optional): Area within which connections are made. Defaults to None.
        max_length (float, optional): Max edge length. Defaults to None.
    """
    if within is None:
        # If not provided, create a box from the maximum bounds
        xmin = min(
            network._mesh1d.mesh1d_node_x.min(), network._mesh2d.mesh2d_node_x.min()
        )
        xmax = max(
            network._mesh1d.mesh1d_node_x.max(), network._mesh2d.mesh2d_node_x.max()
        )
        ymin = min(
            network._mesh1d.mesh1d_node_y.min(), network._mesh2d.mesh2d_node_y.min()
        )
        ymax = max(
            network._mesh1d.mesh1d_node_y.max(), network._mesh2d.mesh2d_node_y.max()
        )

        within = box(xmin, ymin, xmax, ymax)

    # If a 'within' polygon was provided, convert it to a geometrylist
    geometrylist = _geomlist_from_polygon(within)

    # Get the nodes for the specific branch ids
    node_mask = network._mesh1d.get_node_mask(branchids)

    # Get the already present links. These are not filtered on length
    npresent = len(network._link1d2d.link1d2d)

    # Generate links
    network._link1d2d._link_from_1d_to_2d(node_mask, polygon=geometrylist)

    # Filter the links that are longer than the max distance
    id1d = network._link1d2d.link1d2d[npresent:, 0]
    id2d = network._link1d2d.link1d2d[npresent:, 1]
    nodes1d = np.stack(
        [network._mesh1d.mesh1d_node_x[id1d], network._mesh1d.mesh1d_node_y[id1d]],
        axis=1,
    )
    faces2d = np.stack(
        [network._mesh2d.mesh2d_face_x[id2d], network._mesh2d.mesh2d_face_y[id2d]],
        axis=1,
    )
    lengths = np.hypot(nodes1d[:, 0] - faces2d[:, 0], nodes1d[:, 1] - faces2d[:, 1])
    keep = np.concatenate(
        [np.arange(npresent), np.nonzero(lengths < max_length)[0] + npresent]
    )
    _filter_links_on_idx(network, keep)

links1d2d_add_links_2d_to_1d_embedded(network: Network, branchids: list[str] = None, within: Polygon | MultiPolygon = None) -> None

Generates links from 2d to 1d, where the 2d mesh intersects the 1d mesh: the 'embedded' links.

To find the intersecting cells in an efficient way, we follow we the next steps. 1) Get the maximum length of a face edge. 2) Buffer the branches with this length. 3) Find all face nodes within this buffered geometry. 4) Check for each of the corresponding faces if it crossed the branches.

Parameters:

Name	Type	Description	Default
network	Network	Network in which the links are made. Should contain a 1d and 2d mesh	required
branchids	List[str]	List is branch id's for which the connections are made. Defaults to None.	None
within	Union[Polygon, MultiPolygon]	Clipping polygon for 2d mesh that is. Defaults to None.	None

Source code in hydrolib/dhydamo/geometry/mesh.py

def links1d2d_add_links_2d_to_1d_embedded(
    network: Network,
    branchids: list[str] = None,
    within: Polygon | MultiPolygon = None,
) -> None:
    """Generates links from 2d to 1d, where the 2d mesh intersects the 1d mesh: the 'embedded' links.

    To find the intersecting cells in an efficient way, we follow we the next steps. 1) Get the
    maximum length of a face edge. 2) Buffer the branches with this length. 3) Find all face nodes
    within this buffered geometry. 4) Check for each of the corresponding faces if it crossed the
    branches.

    Args:
        network (Network): Network in which the links are made. Should contain a 1d and 2d mesh
        branchids (List[str], optional): List is branch id's for which the connections are made. Defaults to None.
        within (Union[Polygon, MultiPolygon], optional): Clipping polygon for 2d mesh that is. Defaults to None.

    """
    # Get the max edge distance
    nodes2d = np.stack(
        [network._mesh2d.mesh2d_node_x, network._mesh2d.mesh2d_node_y], axis=1
    )
    edge_node_crds = nodes2d[network._mesh2d.mesh2d_edge_nodes]

    diff = edge_node_crds[:, 0, :] - edge_node_crds[:, 1, :]
    maxdiff = np.hypot(diff[:, 0], diff[:, 1]).max()

    # Create multilinestring from branches
    nodes1d = np.stack(
        [network._mesh1d.mesh1d_node_x, network._mesh1d.mesh1d_node_y], axis=1
    )

    # Create a prepared multilinestring of the 1d network, to check for intersections
    mls = MultiLineString(nodes1d[network._mesh1d.mesh1d_edge_nodes].tolist())
    mls_prep = prep(mls)

    # Buffer the branches with the max cell distances
    area = mls.buffer(maxdiff)

    # If a within polygon is provided, clip the buffered area with this polygon.
    if within is not None:
        area = area.intersection(within)

    # Create an array with 2d facecenters and check which intersect the (clipped) area
    faces2d = np.stack(
        [network._mesh2d.mesh2d_face_x, network._mesh2d.mesh2d_face_y], axis=1
    )    
    mpgl = mk.GeometryList(*faces2d.T.copy())
    idx = np.zeros(len(faces2d), dtype=bool)

    for subarea in common.as_polygon_list(area):
        subarea = GeometryList.from_polygon(subarea)
        sub_idx = network.meshkernel.polygon_get_included_points(subarea, mpgl)
        idx |= sub_idx.values == 1

    # Check for each of the remaining faces, if it actually crosses the branches
    nodes2d = np.stack(
        [network._mesh2d.mesh2d_node_x, network._mesh2d.mesh2d_node_y], axis=1
    )
    where = np.nonzero(idx)[0]
    for i, face_idxs in enumerate(network._mesh2d.mesh2d_face_nodes[idx]):
        # nodes2d always has 4 columns, but a triangular cell only has 3 nodes.
        # In that case, one of the indices will be invalid.
        face_crds = nodes2d[face_idxs[face_idxs >= 0]]

        # Check if the face crosses a branch
        if not mls_prep.intersects(LineString(face_crds)):
            idx[where[i]] = False

    # Use the remaining points to create the links
    multipoint = mk.GeometryList(
        x_coordinates=faces2d[idx, 0], y_coordinates=faces2d[idx, 1]
    )

    # Get the nodes for the specific branch ids
    node_mask = network._mesh1d.get_node_mask(branchids)

    # Generate links
    network._link1d2d._link_from_2d_to_1d_embedded(node_mask, polygons=multipoint)

links1d2d_add_links_2d_to_1d_lateral(network: Network, dist_factor: float | None = 2.0, branchids: list[str] = None, within: Polygon | MultiPolygon = None, max_length: float = np.inf) -> None

Generate 1d2d links from the 2d mesh to the 1d mesh, with a lateral connection. If a link is kept, is determined based on the distance between the face center and the intersection with the 2d mesh exterior. By default, links with an intersection distance larger than 2 times the center to edge distance of the cell, are removed. Note that for a square cell with a direct link out of the cell (without passing any other cells) this max distance is sqrt(2) = 1.414. The default value of 2 provides some flexibility. Note that a link with more than 1 intersection with the 2d mesh boundary is removed anyway.

Furthermore: - Branch ids can be specified to connect only specific branches. - A 'within' polygon can be given to only connect 2d cells within this polygon. - A max link length can be given to limit the link length.

Parameters:

Name	Type	Description	Default
network	Network	Network in which the links are made. Should contain a 1d and 2d mesh	required
dist_factor	Union[float, None]	Factor to determine which links are kept (see description above). Defaults to 2.0.	2.0
branchids	List[str]	List is branch id's for which the conncetions are made. Defaults to None.	None
within	Union[Polygon, MultiPolygon]	Clipping polygon for 2d mesh that is. Defaults to None.	None
max_length	float	Max edge length. Defaults to None.	inf

Source code in hydrolib/dhydamo/geometry/mesh.py

def links1d2d_add_links_2d_to_1d_lateral(
    network: Network,
    dist_factor: float | None = 2.0,
    branchids: list[str] = None,
    within: Polygon | MultiPolygon = None,
    max_length: float = np.inf,
) -> None:
    """Generate 1d2d links from the 2d mesh to the 1d mesh, with a lateral connection.
    If a link is kept, is determined based on the distance between the face center and
    the intersection with the 2d mesh exterior. By default, links with an intersection
    distance larger than 2 times the center to edge distance of the cell, are removed.
    Note that for a square cell with a direct link out of the cell (without passing any
    other cells) this max distance is sqrt(2) = 1.414. The default value of 2 provides
    some flexibility. Note that a link with more than 1 intersection with the 2d mesh
    boundary is removed anyway.

    Furthermore:
    - Branch ids can be specified to connect only specific branches.
    - A 'within' polygon can be given to only connect 2d cells within this polygon.
    - A max link length can be given to limit the link length.

    Args:
        network (Network): Network in which the links are made. Should contain a 1d and 2d mesh
        dist_factor (Union[float, None], optional): Factor to determine which links are kept (see description above). Defaults to 2.0.
        branchids (List[str], optional): List is branch id's for which the conncetions are made. Defaults to None.
        within (Union[Polygon, MultiPolygon], optional): Clipping polygon for 2d mesh that is. Defaults to None.
        max_length (float, optional): Max edge length. Defaults to None.
    """

    geometrylist = network.meshkernel.mesh2d_get_mesh_boundaries_as_polygons()
    mpboundaries = GeometryList(**geometrylist.__dict__).to_geometry()
    if within is not None:
        # If a 'within' polygon was provided, get the intersection with the meshboundaries
        # and convert it to a geometrylist
        # Note that the provided meshboundaries is a (list of) polygon(s). Holes are provided
        # as polygons as well, which dont make it a valid MultiPolygon
        if isinstance(mpboundaries, Polygon):
            geom = MultiPolygon([mpboundaries.intersection(within)])
            geometrylist = GeometryList.from_geometry(geom)
        else:
            geom = [geom.intersection(within) for geom in mpboundaries.geoms]
            geom = MultiPolygon(common.as_polygon_list(geom))
            geometrylist = GeometryList.from_geometry(geom)

    # Get the nodes for the specific branch ids
    node_mask = network._mesh1d.get_node_mask(branchids)

    # Get the already present links. These are not filtered subsequently
    present_links = network._link1d2d.link1d2d

    # Generate links
    network._link1d2d._link_from_2d_to_1d_lateral(
        node_mask, polygon=geometrylist, search_radius=max_length
    )

    # If the provided distance factor was None, no further selection is needed, all links are kept.
    if dist_factor is None:
        return

    # Create multilinestring
    if isinstance(mpboundaries, Polygon):
        multilinestring = MultiLineString([mpboundaries.exterior])
    else:
        multilinestring = MultiLineString([poly.exterior for poly in mpboundaries.geoms])

    # Find the links that intersect the boundary close to the origin    
    id1d = network._link1d2d.link1d2d[:, 0]
    id2d = network._link1d2d.link1d2d[:, 1]

    nodes1d = np.stack(
        [network._mesh1d.mesh1d_node_x[id1d], network._mesh1d.mesh1d_node_y[id1d]],
        axis=1,
    )
    faces2d = np.stack(
        [network._mesh2d.mesh2d_face_x[id2d], network._mesh2d.mesh2d_face_y[id2d]],
        axis=1,
    )
    nodes2d = np.stack(
        [network._mesh2d.mesh2d_node_x, network._mesh2d.mesh2d_node_y], axis=1
    )
    nodes2d_idx = network._mesh2d.mesh2d_face_nodes[id2d]
    keep = []
    for i, (node1d, face2d, node2d) in enumerate(
        zip(nodes1d, faces2d, nodes2d_idx)
    ):
        # Triangular grids
        face_node_crds = nodes2d[node2d[node2d > 0.]]                

        # Calculate distance between face edge and face center
        x1 = face_node_crds[:,0]
        y1 = face_node_crds[:,1]
        x0, y0 = face2d[0], face2d[1]
        distance = np.hypot(x1-x0, y1-y0).mean() 

        isect = multilinestring.intersection(LineString([face2d, node1d]))

        # If the intersection is for some reason not a Point of Multipoint, skip it.
        if not isinstance(isect, (Point, MultiPoint)):
            continue

        # Skip the link if it has more than one intersection with the boundary
        isect_list = common.as_point_list(isect)
        if len(isect_list) != 1:
            continue

        # If the distance to the mesh 2d exterior intersection is smaller than
        # the compared distance, keep it.
        dist = np.hypot(*(face2d - isect_list[0].coords[:][0]))
        if dist < dist_factor * distance:
            keep.append(i)

    # Select the remaining links
    _filter_links_on_idx(network, keep, present_links=present_links)

links1d2d_remove_1d_endpoints(network: Network) -> None

Method to remove 1d2d links from end points of the 1d mesh. The GUI will interpret every endpoint as a boundary conditions, which does not allow a 1d 2d link at the same node. To avoid problems with this, use this method.

Source code in hydrolib/dhydamo/geometry/mesh.py

def links1d2d_remove_1d_endpoints(network: Network) -> None:
    """Method to remove 1d2d links from end points of the 1d mesh. The GUI
    will interpret every endpoint as a boundary conditions, which does not
    allow a 1d 2d link at the same node. To avoid problems with this, use
    this method.
    """
    # Can only be done after links have been generated
    if not list(network._mesh1d.network1d_node_id) or not list(
        network._mesh2d.mesh2d_face_nodes
    ):
        logger.warning("1d nodes or 2d faces are not present.")
        return None

    if network._link1d2d.link1d2d.shape[0] == 0:
        logger.warning("No 1d2d-links present.")        
        return None

    # Select mesh1d nodes that are only present in a single edge.
    # link1d2d[:, 0] stores mesh1d node indices, so endpoint detection must
    # use the same mesh1d index space.
    edge_nodes = network._mesh1d.mesh1d_edge_nodes
    edgeid, counts = np.unique(edge_nodes, return_counts=True)
    to_remove = edgeid[counts == 1]

    # Remove links connected to 1d endpoint nodes.
    link_nodes = network._link1d2d.link1d2d[:, 0]
    keep = ~np.isin(link_nodes, to_remove)

    _filter_links_on_idx(network, keep)

links1d2d_remove_within(network: Network, within: Polygon | MultiPolygon) -> None

Remove 1d2d links within a given polygon or multipolygon

Parameters:

Name	Type	Description	Default
network	Network	The network from which the links are removed	required
within	Union[Polygon, MultiPolygon]	The polygon that indicates which to remove	required

Source code in hydrolib/dhydamo/geometry/mesh.py

def links1d2d_remove_within(
    network: Network, within: Polygon | MultiPolygon
) -> None:
    """Remove 1d2d links within a given polygon or multipolygon

    Args:
        network (Network): The network from which the links are removed
        within (Union[Polygon, MultiPolygon]): The polygon that indicates which to remove
    """

    # Create an array with 2d facecenters and 1d nodes, that form the links
    links = network._link1d2d.link1d2d
    nodes1d = np.stack(
        [
            network._mesh1d.mesh1d_node_x[links[:, 0]],
            network._mesh1d.mesh1d_node_y[links[:, 0]],
        ], axis=1
    )
    faces2d = np.stack(
        [
            network._mesh2d.mesh2d_face_x[links[:, 1]],
            network._mesh2d.mesh2d_face_y[links[:, 1]],
        ], axis=1
    )

    # Create GeometryList MultiPoint object (dedupe to avoid MK duplicate-point issues)
    uniq_faces2d, inv_faces2d = np.unique(faces2d, axis=0, return_inverse=True)
    uniq_nodes1d, inv_nodes1d = np.unique(nodes1d, axis=0, return_inverse=True)
    mpgl_faces2d = GeometryList(*uniq_faces2d.T.copy())
    mpgl_nodes1d = GeometryList(*uniq_nodes1d.T.copy())
    idx_faces = np.zeros(len(uniq_faces2d), dtype=bool)
    idx_nodes = np.zeros(len(uniq_nodes1d), dtype=bool)

    # Check which links intersect the provided area
    for polygon in common.as_polygon_list(within):
        subarea = GeometryList.from_geometry(polygon)
        idx_faces |= (
            network.meshkernel.polygon_get_included_points(subarea, mpgl_faces2d).values
            == 1
        )
        idx_nodes |= (
            network.meshkernel.polygon_get_included_points(subarea, mpgl_nodes1d).values
            == 1
        )

    # Remove these links
    idx = idx_faces[inv_faces2d] | idx_nodes[inv_nodes1d]
    keep = ~idx
    _filter_links_on_idx(network, keep)

mesh1d_add_branch_from_linestring(network: Network, linestring: LineString, node_distance: float | int, name: str | None = None, structure_chainage: list[float] | None = None, max_dist_to_struc: float | None = None) -> str

Add branch to 1d mesh, from a LineString geometry. The branch is discretized with the given node distance. The position of a structure can be provided, just like the max distance of a mesh node to a structure

Parameters:

Name	Type	Description	Default
network	Network	Network to which the branch is added	required
linestring	LineString	The geometry of the new branch	required
node_distance	Union[float, int]	Preferred node distance between branch nodes	required
name	Union[str, None]	Name of the branch. If not given, a name is generated. Defaults to None.	None
structure_chainage	Union[List[float], None]	Positions of structures along the branch. Defaults to Union[float, None]=None.	None
max_dist_to_struc	Union[float, None]	Max distance of a mesh point to a structure. Defaults to Union[float, None]=None.	None

Returns:

Name	Type	Description
str	str	name of added branch

Source code in hydrolib/dhydamo/geometry/mesh.py

def mesh1d_add_branch_from_linestring(
    network: Network,
    linestring: LineString,
    node_distance: float | int,
    name: str | None = None,
    structure_chainage: list[float] | None = None,
    max_dist_to_struc: float | None = None,
) -> str:
    """Add branch to 1d mesh, from a LineString geometry.
    The branch is discretized with the given node distance.
    The position of a structure can be provided, just like the max distance
    of a mesh node to a structure

    Args:
        network (Network): Network to which the branch is added
        linestring (LineString): The geometry of the new branch
        node_distance (Union[float, int]): Preferred node distance between branch nodes
        name (Union[str, None], optional): Name of the branch. If not given, a name is generated. Defaults to None.
        structure_chainage (Union[List[float], None], optional): Positions of structures along the branch. Defaults to Union[float, None]=None.
        max_dist_to_struc (Union[float, None], optional): Max distance of a mesh point to a structure. Defaults to Union[float, None]=None.

    Returns:
        str: name of added branch
    """

    branch = Branch(geometry=np.array(linestring.coords[:]))
    branch.generate_nodes(
        mesh1d_edge_length=node_distance,
        structure_chainage=structure_chainage,
        max_dist_to_struc=max_dist_to_struc,
    )
    branchid = network.mesh1d_add_branch(branch, name=name)

    return branchid

mesh1d_add_branches_from_gdf(network: Network, branches: gpd.GeoDataFrame, branch_name_col: str, node_distance: float, max_dist_to_struc: float = None, structures=None) -> None

Function to generate branches from geodataframe

Parameters:

Name	Type	Description	Default
network	Network	The network to which the branches are added	required
branches	GeoDataFrame	GeoDataFrame with branches	required
branch_name_col	str	Name of the column in the GeoDataFrame with the branchnames	required
node_distance	float	Preferred 1d mesh distance	required
max_dist_to_struc	float	Maximum distance to structure. Defaults to None.	None
structures	GeoDataFrame	GeoDataFrame with structures. Must contain a column branchid and chainage. Defaults to None.	None

Source code in hydrolib/dhydamo/geometry/mesh.py

def mesh1d_add_branches_from_gdf(
    network: Network,
    branches: gpd.GeoDataFrame,
    branch_name_col: str,
    node_distance: float,
    max_dist_to_struc: float = None,
    structures=None,
) -> None:
    """Function to generate branches from geodataframe

    Args:
        network (Network): The network to which the branches are added
        branches (gpd.GeoDataFrame): GeoDataFrame with branches
        branch_name_col (str): Name of the column in the GeoDataFrame with the branchnames
        node_distance (float): Preferred 1d mesh distance
        max_dist_to_struc (float, optional): Maximum distance to structure. Defaults to None.
        structures (gpd.GeoDataFrame, optional): GeoDataFrame with structures. Must contain a column branchid and chainage. Defaults to None.
    """

    # Create empty dictionary for structure chainage
    structure_chainage = {}

    # If structures are given, collect offsets per branch
    if structures is not None:
        # Get structure data from dfs
        ids_offsets = structures[["branchid", "chainage", "id"]].copy()
        idx = structures["branchid"] != ""
        missing = structures.branchid.isna()
        if missing.any():
            n_missing = int(missing.sum())
            if n_missing > 5:
                logger.warning(
                    "%d structures are not linked to a branch (too many to show).",
                    n_missing,
                )
            else:
                logger.warning(
                    "%d structures (%s) are not linked to a branch.",
                    n_missing,
                    structures.loc[missing, "id"].to_list(),
                )
        ids_offsets = ids_offsets.loc[idx, :]

        # For each branch
        for branchid, group in ids_offsets.groupby("branchid")[["chainage", "id"]]:
            # Check if structures are located at the same offset
            u, c = np.unique(group.chainage, return_counts=True)
            if any(c > 1):
                logger.warning(
                    "Structures %s have the same location.",
                    ", ".join(group.loc[np.isin(group, u[c > 1]), "id"].tolist()),
                )
            # Add to dictionary
            structure_chainage[branchid] = u

    # Loop over all branches, and add structures
    for branchname, geometry in zip(
        branches[branch_name_col].tolist(), branches["geometry"].tolist()
    ):
        # Create branch
        branch = Branch(geometry=np.array(geometry.coords[:]))
        # Generate nodes on branch
        branch.generate_nodes(
            mesh1d_edge_length=node_distance,
            structure_chainage=structure_chainage[branchname]
            if branchname in structure_chainage
            else None,
            max_dist_to_struc=max_dist_to_struc,
        )
        network.mesh1d_add_branch(branch, name=branchname)

mesh1d_order_numbers_from_attribute(branches: gpd.GeoDataFrame, missing: list, order_attribute: str, network: Network, exceptions: list = []) -> list

mesh1d_order_numbers_from_attribute Method to allocate order numbers to branches based on a common attribute.

Parameters:

Name	Type	Description	Default
branches	GeoDataFrame	HyDAMO GeoDataFGame containing branches.	required
missing	list	list of branches that have so far not been allocated a cross section.	required
order_attribute	str	description. Attribute of branches based on which order numbers are allocated.	required
network	Network	FM network geometry to which the order numbers should be written	required
exceptions	list	List of branch code that should not be given an order number. Defaults to None.	[]

Returns:

Name	Type	Description
missing_after_interpolation	list	list of branches that still do not have a cross section

Source code in hydrolib/dhydamo/geometry/mesh.py

def mesh1d_order_numbers_from_attribute(branches: gpd.GeoDataFrame, missing: list, order_attribute:str, network: Network, exceptions:list=[])-> list:
    """
    mesh1d_order_numbers_from_attribute Method to allocate order numbers to branches based on a common attribute.

    Args:
        branches (GeoDataFrame): HyDAMO GeoDataFGame containing branches.
        missing (list): list of branches that have so far not been allocated a cross section.
        order_attribute (str): _description_. Attribute of branches based on which order numbers are allocated.
        network: FM network geometry to which the order numbers should be written
        exceptions (list, optional): List of branch code that should not be given an order number. Defaults to None.

    Returns:
        missing_after_interpolation (list): list of branches that still do not have a cross section
    """
    j = 0
    branches["order"] = np.nan
    for value in branches[order_attribute].unique():
        if value is not None and not all(x in missing for x in branches.loc[branches.loc[:, order_attribute] == value, "code"]):
            branches.loc[branches.loc[:, order_attribute] == value, "order"] = int(j)
            j = j + 1
    for exception in exceptions:
        branches.loc[branches.code == exception, 'order']  = -1

    interpolation = []
    for ordernr in branches.order.unique():
        if ordernr > 0:
            mesh1d_set_branch_order(
                network,
                branches.code[branches.order == ordernr].to_list(),
                idx=int(ordernr),
            )
            interpolation = (
                interpolation + branches.code[branches.order == ordernr].to_list()
            )
    missing_after_interpolation = np.setdiff1d(missing, interpolation)
    return(missing_after_interpolation)

mesh1d_set_branch_order(network: Network, branchids: list, idx: int = None) -> None

Group branch ids so that the cross sections are interpolated along the branch.

Parameters

branchids : list List of branches to group idx : int Order number with which to update a branch

Source code in hydrolib/dhydamo/geometry/mesh.py

def mesh1d_set_branch_order(network: Network, branchids: list, idx: int = None) -> None:
    """
    Group branch ids so that the cross sections are
    interpolated along the branch.

    Parameters
    ----------
    branchids : list
        List of branches to group
    idx : int
        Order number with which to update a branch
    """
    # Get the ids (integers) of the branch names given by the user
    branchidx = np.isin(network._mesh1d.network1d_branch_id, branchids)
    # Get current order
    branchorder = network._mesh1d.network1d_branch_order
    # Update
    if idx is None:
        branchorder[branchidx] = branchorder.max() + 1
    else:
        if not isinstance(idx, int):
            raise TypeError("Expected integer.")
        branchorder[branchidx] = idx
    # Save
    network._mesh1d.network1d_branch_order = branchorder

mesh2d_add_rectilinear(network: Network, polygon: Polygon | MultiPolygon, dx: float, dy: float, deletemeshoption: mk.DeleteMeshOption = mk.DeleteMeshOption.INSIDE_NOT_INTERSECTED) -> None

Add 2d rectilinear mesh to network. A new network is created, clipped, and merged with the existing network.

Parameters:

Name	Type	Description	Default
network	Network	Network object to which the mesh is added	required
polygon	Union[Polygon, MultiPolygon]	Geometry within which the mesh is generated	required
dx	float	Horizontal mesh spacing	required
dy	float	Vertical mesh spacing	required
deletemeshoption	int	Option for clipping mesh. Defaults to 1.	INSIDE_NOT_INTERSECTED

Returns:

Name	Type	Description
_type_	None	description

Source code in hydrolib/dhydamo/geometry/mesh.py

def mesh2d_add_rectilinear(
    network: Network,
    polygon: Polygon | MultiPolygon,
    dx: float,
    dy: float,
    deletemeshoption: mk.DeleteMeshOption=mk.DeleteMeshOption.INSIDE_NOT_INTERSECTED
) -> None:
    """Add 2d rectilinear mesh to network. A new network is created, clipped, and merged
    with the existing network.

    Args:
        network (Network): Network object to which the mesh is added
        polygon (Union[Polygon, MultiPolygon]): Geometry within which the mesh is generated
        dx (float): Horizontal mesh spacing
        dy (float): Vertical mesh spacing
        deletemeshoption (int, optional): Option for clipping mesh. Defaults to 1.

    Returns:
        _type_: _description_
    """

    # Loop over polygons and create a grid for each polygon bounds
    mesh2d_hc = Mesh2d(meshkernel=mk.MeshKernel())
    polygon_list = GeometryList.from_geometry(polygon)
    for part in polygon_list.geoms:
        xmin, ymin, xmax, ymax = part.bounds
        rows = int((ymax - ymin) / dy)
        columns = int((xmax - xmin) / dx)
        if rows > 0 and columns > 0:
            mesh2d_hc.create_rectilinear(extent=part.bounds, dx=dx, dy=dy)

    # Store present 2d mesh (to be able to add)
    existing_mesh2d = network._mesh2d.get_mesh2d()

    # Create new network
    mesh2d_output = mesh2d_hc.get_mesh2d()
    network._mesh2d._set_mesh2d(
        mesh2d_output.node_x,
        mesh2d_output.node_y,
        mesh2d_output.edge_nodes,
    )

    # Clip and clean
    mesh2d_clip(
        network=network,
        polygon=polygon_list,
        deletemeshoption=deletemeshoption,
        inside=False,
    )

    # Merge with existing network
    if existing_mesh2d.node_x.size > 0:
        # Add all variables to existing mesh
        new_mesh2d = network._mesh2d.get_mesh2d()
        network._mesh2d._set_mesh2d(
            np.concatenate([existing_mesh2d.node_x, new_mesh2d.node_x]),
            np.concatenate([existing_mesh2d.node_y, new_mesh2d.node_y]),
            np.concatenate([
                existing_mesh2d.edge_nodes,
                new_mesh2d.edge_nodes + existing_mesh2d.edge_nodes.max() + 1,
            ]),
        )

mesh2d_add_triangular(network: Network, polygon: Polygon | MultiPolygon, edge_length: float = None) -> None

Add triangular mesh to existing network. An orthogonal mesh is generated by the meshkernel, which likely means that the given geometry is not completely filled. The triangle discretization is determined based on the coordinates on the boundary of the provided geometry. Giving an edge_length will discretize the polygon for you, but you can also do this yourself.

Parameters:

Name	Type	Description	Default
network	Network	Network object to which the mesh is added	required
polygon	Union[Polygon, MultiPolygon]	Geometry within which the mesh is generated	required
edge_length	float	Distance for which the polygon boundary is discretized (by approximation). Defaults to None.	None

Source code in hydrolib/dhydamo/geometry/mesh.py

def mesh2d_add_triangular(
    network: Network, polygon: Polygon | MultiPolygon, edge_length: float = None
) -> None:
    """Add triangular mesh to existing network. An orthogonal mesh is generated by the
    meshkernel, which likely means that the given geometry is not completely filled. The
    triangle discretization is determined based on the coordinates on the boundary of the
    provided geometry. Giving an edge_length will discretize the polygon for you, but
    you can also do this yourself.

    Args:
        network (Network): Network object to which the mesh is added
        polygon (Union[Polygon, MultiPolygon]): Geometry within which the mesh is generated
        edge_length (float, optional): Distance for which the polygon boundary is discretized (by approximation). Defaults to None.
    """

    meshkernel = network._mesh2d.meshkernel
    for polygon in common.as_polygon_list(polygon):
        # Interpolate coordinates on polygon with edge_length distance
        if edge_length is not None:
            polygon = common.interp_polygon(polygon, dist=edge_length)

        # Add triangular mesh within polygon
        meshkernel.mesh2d_make_triangular_mesh_from_polygon(
            _geomlist_from_polygon(polygon),
            scale_factor=1,
            )

mesh2d_altitude_from_raster(network, rasterpath, where: RasterStatPosition = 'face', stat='mean', fill_option: FillOption = 'fill_value', fill_value=None)

Method to determine level of nodes

This function works faster for large amounts of cells, since it does not draw polygons but checks for nearest neighbours (voronoi) based on interpolation.

Note that the raster is not clipped. Any values outside the bounds are also taken into account.

Source code in hydrolib/dhydamo/geometry/mesh.py

def mesh2d_altitude_from_raster(
    network,
    rasterpath,
    where: RasterStatPosition = "face",
    stat="mean",
    fill_option: FillOption = "fill_value",
    fill_value=None,
):
    """
    Method to determine level of nodes

    This function works faster for large amounts of cells, since it does not
    draw polygons but checks for nearest neighbours (voronoi) based
    on interpolation.

    Note that the raster is not clipped. Any values outside the bounds are
    also taken into account.
    """

    if isinstance(fill_option, str):
        fill_option = FillOption(fill_option)

    if isinstance(where, str):
        where = RasterStatPosition(where)

    # Select points on faces or nodes
    logger.info("Creating GeoDataFrame of cell faces.")

    # Get coordinates where z-values are derived or centered
    xy = np.stack(
        [
            getattr(network._mesh2d, f"mesh2d_{where.value}_x"),
            getattr(network._mesh2d, f"mesh2d_{where.value}_y"),
        ],
        axis=1,
    )

    # Create cells as polygons
    xy_facenodes = np.stack(
        [network._mesh2d.mesh2d_node_x, network._mesh2d.mesh2d_node_y], axis=1
    )
    cells = network._mesh2d.mesh2d_face_nodes
    nodatavalue = np.iinfo(cells.dtype).min
    indices = cells != nodatavalue
    indices = cells != -2147483648
    cells = [xy_facenodes[cell[index]] for cell, index in zip(cells, indices)]
    facedata = gpd.GeoDataFrame(geometry=[Polygon(cell) for cell in cells])

    if where == RasterStatPosition.FACE:
        # Get raster statistics
        facedata.index = np.arange(len(xy), dtype=np.uint32) + 1
        facedata["crds"] = [cell for cell in cells]

        df = rasterstats.raster_stats_fine_cells(rasterpath, facedata, stats=[stat])
        # Get z values
        zvalues = df[stat].values

    elif where == RasterStatPosition.NODE:
        logger.info(
            "Generating voronoi polygons around cell centers for determining raster statistics."
        )

        facedata = spatial.get_voronoi_around_nodes(xy, facedata)
        # Get raster statistics
        df = rasterstats.raster_stats_fine_cells(rasterpath, facedata, stats=[stat])
        # Get z values
        zvalues = df[stat].values

    isnan = np.isnan(zvalues)
    if isnan.any():
        # If interpolation or fill_option, but all are NaN, raise error.
        if isnan.all() and fill_option in [FillOption.NEAREST, FillOption.INTERPOLATE]:
            raise ValueError(
                "Only NaN values found, interpolation or nearest not possible."
            )

        # Fill fill_option values
        if fill_option == FillOption.FILL_VALUE:
            if fill_value is None:
                raise ValueError("Provide a fill_value (keyword argument).")
            # With default value
            zvalues[isnan] = fill_value

        elif fill_option == FillOption.NEAREST:
            # By looking for the nearest value in the grid
            # Create a KDTree of the known points
            tree = KDTree(data=xy[~isnan])
            idx = tree.query(x=xy[isnan])[1]
            zvalues[isnan] = zvalues[~isnan][idx]

        elif fill_option == FillOption.INTERPOLATE:
            if fill_value is None:
                raise ValueError(
                    "Provide a fill_value (keyword argument) to fill values that cannot be interpolated."
                )
            # By interpolating
            isnan = np.isnan(zvalues)
            interp = LinearNDInterpolator(
                xy[~isnan], zvalues[~isnan], fill_value=fill_value
            )
            zvalues[isnan] = interp(xy[isnan])

    # Set values to mesh geometry
    setattr(network._mesh2d, f"mesh2d_{where.value}_z", zvalues)

mesh2d_clip(network: Network, polygon: GeometryList | (Polygon | MultiPolygon), deletemeshoption: mk.DeleteMeshOption = mk.DeleteMeshOption.INSIDE_AND_INTERSECTED, inside=True) -> None

Clip the mesh (currently implemented for 2d) and clean remaining hanging edges.

Parameters:

Name	Type	Description	Default
network	Network	Network for which the mesh is clipped	required
polygon	Union[GeometryList, Union[Polygon, MultiPolygon]]	Polygon within which the mesh is clipped	required
deletemeshoption	int	Options for deleting nodes inside/outside polygon. Defaults to 1.	INSIDE_AND_INTERSECTED
inside	bool	Whether to clip inside or outside the polygon. Defaults to True.	True

Source code in hydrolib/dhydamo/geometry/mesh.py

def mesh2d_clip(
    network: Network,
    polygon: GeometryList | (Polygon | MultiPolygon),
    deletemeshoption: mk.DeleteMeshOption = mk.DeleteMeshOption.INSIDE_AND_INTERSECTED,
    inside=True,
) -> None:
    """Clip the mesh (currently implemented for 2d) and clean remaining hanging edges.

    Args:
        network (Network): Network for which the mesh is clipped
        polygon (Union[GeometryList, Union[Polygon, MultiPolygon]]): Polygon within which the mesh is clipped
        deletemeshoption (int, optional): Options for deleting nodes inside/outside polygon. Defaults to 1.
        inside (bool, optional): Whether to clip inside or outside the polygon. Defaults to True.
    """

    if isinstance(polygon, GeometryList):
        geo = polygon.to_geometry()
        if not isinstance(geo, (Polygon, MultiPolygon)):
            raise TypeError(
                f"Expected to provided geometrylist to be interpreted as Polygon or MultiPolygon, not a {type(geo)}."
            )
    elif isinstance(polygon, (Polygon, MultiPolygon)):
        polygon = GeometryList.from_geometry(polygon)

    if inside:
        # Clip each single polygon
        for part in polygon.geoms:
            network.mesh2d_clip_mesh(_geomlist_from_polygon(part), deletemeshoption, inside)
    else:
        # Keep exterior
        clip_pol_ext = _geomlist_from_multipolygon(polygon, "exterior")
        if clip_pol_ext.x_coordinates.size > 0:
            network._mesh2d.meshkernel.mesh2d_delete(
                geometry_list=clip_pol_ext,
                delete_option=deletemeshoption,
                invert_deletion=True,
            )
        # # Clip interior
        clip_pol_int = _geomlist_from_multipolygon(polygon, "interior")
        if clip_pol_int.x_coordinates.size > 0:
            network._mesh2d.meshkernel.mesh2d_delete(
                geometry_list=clip_pol_int,
                delete_option=deletemeshoption,
                invert_deletion=False,
            )

    # Remove hanging edges
    network._mesh2d.meshkernel.mesh2d_delete_hanging_edges()    
    network._mesh2d.meshkernel.mesh2d_get()

mesh2d_refine(network: Network, polygon: Polygon | MultiPolygon, steps: int, refine_parameters: dict = None) -> None

Refine mesh 2d within (list of) polygon or multipolygon, with a certain number of refinement steps.

Parameters:

Name	Type	Description	Default
network	Network	Network for which the mesh is clipped	required
polygon	Union[GeometryList, Union[Polygon, MultiPolygon]]	Polygon within which the mesh is clipped	required
steps	int	Number of steps in the refinement	required
parameters	MeshRefinementParameters	parameters to be used in the refinement. Defaults to None.	required

Source code in hydrolib/dhydamo/geometry/mesh.py

def mesh2d_refine(
    network: Network,    
    polygon: Polygon | MultiPolygon,
    steps: int,
    refine_parameters: dict=None,
) -> None:
    """Refine mesh 2d within (list of) polygon or multipolygon, with a certain
    number of refinement steps.

    Args:
        network (Network): Network for which the mesh is clipped
        polygon (Union[GeometryList, Union[Polygon, MultiPolygon]]): Polygon within which the mesh is clipped
        steps (int): Number of steps in the refinement
        parameters (MeshRefinementParameters, optional): parameters to be used in the refinement. Defaults to None. 
    """
    parameters = mk.MeshRefinementParameters(   
                refine_intersected=True,
                use_mass_center_when_refining=True,
                min_edge_size=1.0,
                refinement_type=2,
                connect_hanging_nodes=True,
                account_for_samples_outside_face=True,
                max_refinement_iterations=1,
                smoothing_iterations=5,
                max_courant_time=120.0,
                directional_refinement=False
    )
    if refine_parameters is not None:
        for key in refine_parameters.keys():
            setattr(parameters, key, refine_parameters[key])

    for polygon in common.as_polygon_list(polygon):
        network.mesh2d_refine_mesh(
            _geomlist_from_polygon(polygon),
            level=steps,            
            parameters=parameters)                      

mesh2d_gridgeom

Mesh2D_GG

Source code in hydrolib/dhydamo/geometry/mesh2d_gridgeom.py

class Mesh2D_GG:

    def __init__(self, dll_path=None):
        self.fill_value_z = -999.0
        self.missing_z_value = None
        self.dll_path = dll_path
        self.meshgeomdim = meshgeomdim(pointer(c_char()), 2, 0, 0, 0, 0, -1, -1, -1, -1, -1, 0)
        self.meshgeom = meshgeom(self.meshgeomdim)

    def clip_nodes(self, xnodes, ynodes, edge_nodes, polygon, keep_outside=False):
        """
        Method that returns a set of nodes that is clipped within a given polygon.

        Parameters
        ----------
        xnodes : 
            X-coordinates of mesh nodes
        ynodes : 
            Y-coordinates of mesh nodes
        edge_nodes : 
            node ids of edges. These are assumed in 1-based numbering
        polygon : Polygon, Multipolygon or list of

        """

        # Find nodes within one of the polygons
        logger.info('Selecting nodes within polygon.')
        index = np.zeros(len(xnodes), dtype=bool)
        for poly in geometry.as_polygon_list(polygon):
            index = index | geometry.points_in_polygon(np.c_[xnodes, ynodes], poly)
        if keep_outside:
            index = ~index

        # Filter edge nodes, keep the ones that are within a polygon with both points
        nodeids_in_polygon = np.where(index)[0] + 1
        edge_nodes = edge_nodes[np.isin(edge_nodes, nodeids_in_polygon).all(axis=1)]

        # Remove edges that intersect the clipping polygons. This is needed for when the clipping polygon is narrow,
        # So it crosses an edge without containing one of the edge nodes
        logger.info('Creating LineString for each edge.')
        edges = [LineString(line) for line in np.c_[xnodes, ynodes][edge_nodes - 1]]

        logger.info('Checking for intersections.')
        isect = np.zeros(len(edge_nodes), dtype=bool)
        for poly in geometry.as_polygon_list(polygon):
            prepped = prep(poly.exterior)
            isect = isect | np.array([prepped.intersects(e) for e in edges])
            for interior in poly.interiors:
                prepped = prep(interior)
                isect = isect | np.array([prepped.intersects(e) for e in edges])
        edge_nodes = edge_nodes[~isect]

        # Remove edges that have a 1 count, until none are left
        logger.info('Remove edges with only a single node within the clip area.')
        unique, count = np.unique(edge_nodes, return_counts=True)
        while any(count == 1):
            count1nodes = unique[count == 1]
            edge_nodes = edge_nodes[~np.isin(edge_nodes, count1nodes).any(axis=1)]
            unique, count = np.unique(edge_nodes, return_counts=True)

        # Select remaining nodes
        xnodes = xnodes[unique - 1]
        ynodes = ynodes[unique - 1]

        # Make mapping for new id's
        new_id_mapping = {old_id: new_id + 1 for new_id, old_id in enumerate(unique)}
        edge_nodes = np.reshape([new_id_mapping[old_id] for old_id in edge_nodes.ravel()], edge_nodes.shape)

        return xnodes, ynodes, edge_nodes

    def clip_mesh_by_polygon(self, polygon, outside=True):
        """
        Removes part of the existing mesh within the given polygon.

        Parameters
        ----------
        polygon : list, Polygon, MultiPolygon
            A polygon, multi-polygon or list of multiple polygons or multipolygons.
            The faces within the polygons are removed.
        """
        logger.info('Clipping 2D mesh by polygon.')

        # Get xnodes, ynodes and edge_nodes
        xnodes = self.meshgeom.get_values('nodex', as_array=True)
        ynodes = self.meshgeom.get_values('nodey', as_array=True)
        edge_nodes = self.meshgeom.get_values('edge_nodes', as_array=True)

        # Clip
        logger.info('Clipping nodes.')
        xnodes, ynodes, edge_nodes = self.clip_nodes(xnodes, ynodes, edge_nodes, polygon, keep_outside=outside)

        # Update dimensions
        self.meshgeomdim.numnode = len(xnodes)
        self.meshgeomdim.numedge = len(edge_nodes)

        # Add nodes and links
        self.meshgeom.set_values('nodex', xnodes)
        self.meshgeom.set_values('nodey', ynodes)
        self.meshgeom.set_values('edge_nodes', np.ravel(edge_nodes).tolist())

        # Determine what the cells are
        logger.info('Finding cell faces within clipped nodes.')
        self._find_cells(self.meshgeom, dll_path=self.dll_path)

        # Clean nodes, this function deletes nodes based on no longer existing cells
        face_nodes = self.meshgeom.get_values('face_nodes', as_array=True)
        logger.info('Removing incomplete faces.')
        cleaned = self.clean_nodes(xnodes, ynodes, edge_nodes, face_nodes)

        # If cleaning leaves no whole cells, return None
        if cleaned is None:
            raise ValueError('Clipping the grid does not leave any nodes.')
        # Else, get the arguments from the returned tuple
        else:
            xnodes, ynodes, edge_nodes, face_nodes = cleaned

        logger.info('Update mesh with clipped nodes, edges and faces.')
        # Update dimensions
        self.meshgeomdim.numnode = len(xnodes)
        self.meshgeomdim.numedge = len(edge_nodes)

        # Add nodes and links
        self.meshgeom.set_values('nodex', xnodes)
        self.meshgeom.set_values('nodey', ynodes)
        self.meshgeom.set_values('edge_nodes', np.ravel(edge_nodes).tolist())
        self.meshgeom.set_values('face_nodes', np.ravel(face_nodes).tolist())

    def clean_nodes(self, xnodes, ynodes, edge_nodes, face_nodes):
        """
        Method to clean the nodes. Edges that do not form a cell are deleted.
        Alternative method based on checking tuple pairs
        """

        # Select nodes the occur in any face
        node_selection = np.unique(face_nodes)
        node_selection = node_selection[node_selection!=-999]
        # If no nodes (because no faces), return None
        if not node_selection.any():
            return None

        # Remove all edges that are not in this selection
        edge_nodes = edge_nodes[np.isin(edge_nodes, node_selection).all(axis=1)]

        # Remove all segments that are not part of any face
        # Direction is irrelevant, sort to compare in ascending order
        # Convert to a contiguous array for faster comparison
        sorted_edges = np.ascontiguousarray(np.sort(edge_nodes, axis=1)).view(np.dtype((np.void, edge_nodes.dtype.itemsize * edge_nodes.shape[1])))
        idx = np.zeros(len(sorted_edges), dtype=bool)

        # Get combinations of face nodes that can form an edge
        combs = list(combinations(range(len(face_nodes[0])), 2))
        for comb in combs:
            # Get the first, second, third, etc. possible edge of every face
            arr = face_nodes[:, comb]
            # Remove nodata (so for example the non-existing fourth edge of a triangular face)
            arr = arr[(arr != -999).all(axis=1)]
            arr = np.sort(arr, axis=1)
            # Convert to a contiguous array for faster comparison
            sorted_face_edges = np.ascontiguousarray(arr).view(np.dtype((np.void, arr.dtype.itemsize * arr.shape[1])))
            # Check if face edges are in edges
            idx |= np.isin(sorted_edges, sorted_face_edges)[:, 0]

        # Select the nodes that are present in faces
        edge_nodes = edge_nodes[idx]

        xnodes = xnodes[node_selection - 1]
        ynodes = ynodes[node_selection - 1]

            # Make mapping for new id's
        new_id_mapping = {old_id: new_id + 1 for new_id, old_id in enumerate(node_selection)}    

        edge_nodes = np.reshape([-999 if old_id == -999 else new_id_mapping[old_id] for old_id in edge_nodes.ravel()], edge_nodes.shape)
        face_nodes = np.reshape([-999 if old_id == -999 else new_id_mapping[old_id] for old_id in face_nodes.ravel()], face_nodes.shape)

        return xnodes, ynodes, edge_nodes, face_nodes

    @staticmethod
    def _find_cells(geometries, maxnumfacenodes=None, dll_path=None):
        """
        Determine what the cells are in a grid.
        """

        dimensions = geometries.meshgeomdim
        os.add_dll_directory(dll_path)
        wrapperGridgeom = CDLL(os.path.join(os.path.dirname(__file__), '..', 'resources', 'lib', 'gridgeom', 'gridgeom.dll'))
        ierr = wrapperGridgeom.ggeo_deallocate()
        assert ierr == 0

        start_index = c_int(1)

        meshdimout = meshgeomdim()
        meshout = meshgeom(meshdimout)

        ierr = wrapperGridgeom.ggeo_find_cells(
            byref(dimensions),
            byref(geometries),
            byref(meshdimout),
            byref(meshout),
            byref(start_index)
        )
        assert ierr == 0

        dimensions.numface = meshdimout.numface
        if maxnumfacenodes is None:
            dimensions.maxnumfacenodes = meshdimout.maxnumfacenodes
        else:
            dimensions.maxnumfacenodes = maxnumfacenodes

        # Allocate
        for mesh in [geometries, meshout]:
            for var in ['facex', 'facey', 'face_nodes']:
                mesh.allocate(var)

        ierr = wrapperGridgeom.ggeo_find_cells(
            byref(dimensions),
            byref(geometries),
            byref(meshdimout),
            byref(meshout),
            byref(start_index)
        )
        assert ierr == 0

        # Copy content to self meshgeom
        for var in ['facex', 'facey', 'face_nodes']:
            # Create array of allocated size. In case of deviating maxnumfacenodes
            if (maxnumfacenodes is not None) and (var == 'face_nodes'):
                # Create an empty integer. Empty values should only be present in the integer arrays
                arr = np.full(geometries.get_dimensions(var), -999, dtype=int)
                # Get array
                arr[:, :meshout.meshgeomdim.maxnumfacenodes] = meshout.get_values(var, as_array=True)
                # Set array
                geometries.set_values(var, arr.ravel())
            else:
                # Set array
                geometries.set_values(var, meshout.get_values(var))


        ierr = wrapperGridgeom.ggeo_deallocate()

    def altitude_constant(self, constant, where='face'):

        zvalues = np.ones(getattr(self.meshgeomdim, f'num{where}')) * constant
        self.meshgeom.set_values(f'{where}z', zvalues)

    def altitude_from_raster(self, rasterpath, where='face', stat='mean', missing='default'):
        """
        Method to determine level of nodes

        This function works faster for large amounts of cells, since it does not
        draw polygons but checks for nearest neighbours (voronoi) based
        on interpolation.

        Note that the raster is not clipped. Any values outside the bounds are
        also taken into account.
        """

        # Select points on faces or nodes
        logger.info('Creating GeoDataFrame of cell faces.')
        xy = np.c_[self.meshgeom.get_values(f'{where}x'), self.meshgeom.get_values(f'{where}y')]
        cells = self.meshgeom.get_faces()
        facedata = gpd.GeoDataFrame(geometry=[Polygon(cell) for cell in cells])   

        if where == 'face':
            # Get raster statistics
            facedata.index = np.arange(len(xy), dtype=np.uint32) + 1
            # facedata.to_file('test_face.shp')
            facedata['crds'] = [cell for cell in cells]

            df = geometry.raster_stats_fine_cells(rasterpath, facedata, stats=[stat])
            # Get z values
            zvalues = df[stat].values

        elif where == 'node':

            logger.info('Generating voronoi polygons around cell centers for determining raster statistics.')

            # Creat voronoi polygon
            # Add border to limit polygons
            border = box(xy[:, 0].min(), xy[:, 1].min(), xy[:, 0].max(), xy[:, 1].max()).buffer(1000).exterior
            borderpts = [border.interpolate(dist).coords[0] for dist in np.linspace(0, border.length, max(20, border.length//100))]
            vor = Voronoi(points=xy.tolist()+borderpts)
            clippoly = facedata.union_all()
            # Get lines
            lines = []
            for poly in geometry.as_polygon_list(clippoly):
                lines.append(poly.exterior)
                lines.extend([line for line in poly.interiors])
            linesprep = prep(MultiLineString(lines))
            clipprep = prep(clippoly)

            # Collect polygons
            data = []
            for (pr, (i, pt)) in zip(vor.point_region, enumerate(xy)):
                region = vor.regions[pr]
                if pr == -1:
                    break
                while -1 in region:
                    region.remove(-1)
                if len(region) < 3:
                    continue
                crds = vor.vertices[region]
                if clipprep.intersects(Point(pt)):
                    poly = Polygon(crds)
                    if linesprep.intersects(poly):
                        poly = poly.intersection(clippoly)
                        if isinstance(poly, MultiPolygon):
                            poly = poly.buffer(0.001)
                        if isinstance(poly, MultiPolygon):
                            logger.warning('Got multipolygon when clipping voronoi polygon. Only adding coordinates for largest of the polygons.')
                            poly = poly[np.argmax([p.area for p in geometry.as_polygon_list(poly)])]
                        crds = np.vstack(poly.exterior.coords[:])
                    data.append({'geometry': poly, 'crds': crds})

            # Limit to model extend
            facedata = gpd.GeoDataFrame(data)
            facedata.index=np.arange(len(xy), dtype=np.uint32) + 1
            # facedata.drop('crds', axis=1).to_file('test_node.shp')
            # Get raster statistics
            df = geometry.raster_stats_fine_cells(rasterpath, facedata, stats=[stat])
            # Get z values
            zvalues = df[stat].values

        else:
            raise NotImplementedError()


        # If there are no NaN's, return the answer
        isnan = np.isnan(zvalues)
        if not isnan.any():
            # Set values to mesh geometry
            self.meshgeom.set_values(f'{where}z', zvalues)
            return None

        # If interpolation or missing, but all are NaN, raise error.
        elif isnan.all() and missing in ['nearest', 'interpolation']:
            raise ValueError('Only NaN values found, interpolation or nearest not possible.')

        # Fill missing values
        if missing.lower() == 'default':
            # With default value
            zvalues[isnan] = self.fill_value_z

        elif isinstance(missing, (float, int)):
            # With a given number
            zvalues[isnan] = missing

        elif missing.lower() == 'nearest':
            # By looking for the nearest value in the grid
            # Create a KDTree of the known points
            tree = KDTree(data=xy[~isnan])
            idx = tree.query(x=xy[isnan])[1]
            zvalues[isnan] = zvalues[~isnan][idx]

        elif missing.lower() == 'interpolation':
            # By interpolating
            isnan = np.isnan(zvalues)
            interp = LinearNDInterpolator(xy[~isnan], zvalues[~isnan], fill_value=self.fill_value_z)
            zvalues[isnan] = interp(xy[isnan])

        else:
            raise ValueError(f'{missing} not recognized. Choose \'default\', \'nearest\', \'interpolation\' or a number with which to fill the missing data.')

        # Set values to mesh geometry
        self.meshgeom.set_values(f'{where}z', zvalues)

    def set_missing_z_value(self, value):
        self.missing_z_value = value

    def geom_from_netcdf(self, file, only2d=False):
        gridio.from_netcdf_old(self.meshgeom, file, only2d=only2d)
        self._find_cells(self.meshgeom, dll_path=self.dll_path)

    def faces_to_centroid(self):
        """
        The find_cells function does not always return face coordinates
        that are all accepted by the interactor. The coordinates
        are placed on the circumcenters, and placed on the face border
        in case of a coordinate that is outside the cells.

        This function shifts the face coordinates towards the centroid
        (center of gravity) of the cell.
        """
        # Calculate centroids
        centroids = np.vstack([cell.mean(axis=0) for cell in self.meshgeom.get_faces()]).T

        # Set values back to geometry
        self.meshgeom.set_values('facex', centroids[0])
        self.meshgeom.set_values('facey', centroids[1])

altitude_from_raster(rasterpath, where='face', stat='mean', missing='default')

Method to determine level of nodes

This function works faster for large amounts of cells, since it does not draw polygons but checks for nearest neighbours (voronoi) based on interpolation.

Note that the raster is not clipped. Any values outside the bounds are also taken into account.

Source code in hydrolib/dhydamo/geometry/mesh2d_gridgeom.py

def altitude_from_raster(self, rasterpath, where='face', stat='mean', missing='default'):
    """
    Method to determine level of nodes

    This function works faster for large amounts of cells, since it does not
    draw polygons but checks for nearest neighbours (voronoi) based
    on interpolation.

    Note that the raster is not clipped. Any values outside the bounds are
    also taken into account.
    """

    # Select points on faces or nodes
    logger.info('Creating GeoDataFrame of cell faces.')
    xy = np.c_[self.meshgeom.get_values(f'{where}x'), self.meshgeom.get_values(f'{where}y')]
    cells = self.meshgeom.get_faces()
    facedata = gpd.GeoDataFrame(geometry=[Polygon(cell) for cell in cells])   

    if where == 'face':
        # Get raster statistics
        facedata.index = np.arange(len(xy), dtype=np.uint32) + 1
        # facedata.to_file('test_face.shp')
        facedata['crds'] = [cell for cell in cells]

        df = geometry.raster_stats_fine_cells(rasterpath, facedata, stats=[stat])
        # Get z values
        zvalues = df[stat].values

    elif where == 'node':

        logger.info('Generating voronoi polygons around cell centers for determining raster statistics.')

        # Creat voronoi polygon
        # Add border to limit polygons
        border = box(xy[:, 0].min(), xy[:, 1].min(), xy[:, 0].max(), xy[:, 1].max()).buffer(1000).exterior
        borderpts = [border.interpolate(dist).coords[0] for dist in np.linspace(0, border.length, max(20, border.length//100))]
        vor = Voronoi(points=xy.tolist()+borderpts)
        clippoly = facedata.union_all()
        # Get lines
        lines = []
        for poly in geometry.as_polygon_list(clippoly):
            lines.append(poly.exterior)
            lines.extend([line for line in poly.interiors])
        linesprep = prep(MultiLineString(lines))
        clipprep = prep(clippoly)

        # Collect polygons
        data = []
        for (pr, (i, pt)) in zip(vor.point_region, enumerate(xy)):
            region = vor.regions[pr]
            if pr == -1:
                break
            while -1 in region:
                region.remove(-1)
            if len(region) < 3:
                continue
            crds = vor.vertices[region]
            if clipprep.intersects(Point(pt)):
                poly = Polygon(crds)
                if linesprep.intersects(poly):
                    poly = poly.intersection(clippoly)
                    if isinstance(poly, MultiPolygon):
                        poly = poly.buffer(0.001)
                    if isinstance(poly, MultiPolygon):
                        logger.warning('Got multipolygon when clipping voronoi polygon. Only adding coordinates for largest of the polygons.')
                        poly = poly[np.argmax([p.area for p in geometry.as_polygon_list(poly)])]
                    crds = np.vstack(poly.exterior.coords[:])
                data.append({'geometry': poly, 'crds': crds})

        # Limit to model extend
        facedata = gpd.GeoDataFrame(data)
        facedata.index=np.arange(len(xy), dtype=np.uint32) + 1
        # facedata.drop('crds', axis=1).to_file('test_node.shp')
        # Get raster statistics
        df = geometry.raster_stats_fine_cells(rasterpath, facedata, stats=[stat])
        # Get z values
        zvalues = df[stat].values

    else:
        raise NotImplementedError()


    # If there are no NaN's, return the answer
    isnan = np.isnan(zvalues)
    if not isnan.any():
        # Set values to mesh geometry
        self.meshgeom.set_values(f'{where}z', zvalues)
        return None

    # If interpolation or missing, but all are NaN, raise error.
    elif isnan.all() and missing in ['nearest', 'interpolation']:
        raise ValueError('Only NaN values found, interpolation or nearest not possible.')

    # Fill missing values
    if missing.lower() == 'default':
        # With default value
        zvalues[isnan] = self.fill_value_z

    elif isinstance(missing, (float, int)):
        # With a given number
        zvalues[isnan] = missing

    elif missing.lower() == 'nearest':
        # By looking for the nearest value in the grid
        # Create a KDTree of the known points
        tree = KDTree(data=xy[~isnan])
        idx = tree.query(x=xy[isnan])[1]
        zvalues[isnan] = zvalues[~isnan][idx]

    elif missing.lower() == 'interpolation':
        # By interpolating
        isnan = np.isnan(zvalues)
        interp = LinearNDInterpolator(xy[~isnan], zvalues[~isnan], fill_value=self.fill_value_z)
        zvalues[isnan] = interp(xy[isnan])

    else:
        raise ValueError(f'{missing} not recognized. Choose \'default\', \'nearest\', \'interpolation\' or a number with which to fill the missing data.')

    # Set values to mesh geometry
    self.meshgeom.set_values(f'{where}z', zvalues)

clean_nodes(xnodes, ynodes, edge_nodes, face_nodes)

Method to clean the nodes. Edges that do not form a cell are deleted. Alternative method based on checking tuple pairs

Source code in hydrolib/dhydamo/geometry/mesh2d_gridgeom.py

def clean_nodes(self, xnodes, ynodes, edge_nodes, face_nodes):
    """
    Method to clean the nodes. Edges that do not form a cell are deleted.
    Alternative method based on checking tuple pairs
    """

    # Select nodes the occur in any face
    node_selection = np.unique(face_nodes)
    node_selection = node_selection[node_selection!=-999]
    # If no nodes (because no faces), return None
    if not node_selection.any():
        return None

    # Remove all edges that are not in this selection
    edge_nodes = edge_nodes[np.isin(edge_nodes, node_selection).all(axis=1)]

    # Remove all segments that are not part of any face
    # Direction is irrelevant, sort to compare in ascending order
    # Convert to a contiguous array for faster comparison
    sorted_edges = np.ascontiguousarray(np.sort(edge_nodes, axis=1)).view(np.dtype((np.void, edge_nodes.dtype.itemsize * edge_nodes.shape[1])))
    idx = np.zeros(len(sorted_edges), dtype=bool)

    # Get combinations of face nodes that can form an edge
    combs = list(combinations(range(len(face_nodes[0])), 2))
    for comb in combs:
        # Get the first, second, third, etc. possible edge of every face
        arr = face_nodes[:, comb]
        # Remove nodata (so for example the non-existing fourth edge of a triangular face)
        arr = arr[(arr != -999).all(axis=1)]
        arr = np.sort(arr, axis=1)
        # Convert to a contiguous array for faster comparison
        sorted_face_edges = np.ascontiguousarray(arr).view(np.dtype((np.void, arr.dtype.itemsize * arr.shape[1])))
        # Check if face edges are in edges
        idx |= np.isin(sorted_edges, sorted_face_edges)[:, 0]

    # Select the nodes that are present in faces
    edge_nodes = edge_nodes[idx]

    xnodes = xnodes[node_selection - 1]
    ynodes = ynodes[node_selection - 1]

        # Make mapping for new id's
    new_id_mapping = {old_id: new_id + 1 for new_id, old_id in enumerate(node_selection)}    

    edge_nodes = np.reshape([-999 if old_id == -999 else new_id_mapping[old_id] for old_id in edge_nodes.ravel()], edge_nodes.shape)
    face_nodes = np.reshape([-999 if old_id == -999 else new_id_mapping[old_id] for old_id in face_nodes.ravel()], face_nodes.shape)

    return xnodes, ynodes, edge_nodes, face_nodes

clip_mesh_by_polygon(polygon, outside=True)

Removes part of the existing mesh within the given polygon.

Parameters

polygon : list, Polygon, MultiPolygon A polygon, multi-polygon or list of multiple polygons or multipolygons. The faces within the polygons are removed.

Source code in hydrolib/dhydamo/geometry/mesh2d_gridgeom.py

def clip_mesh_by_polygon(self, polygon, outside=True):
    """
    Removes part of the existing mesh within the given polygon.

    Parameters
    ----------
    polygon : list, Polygon, MultiPolygon
        A polygon, multi-polygon or list of multiple polygons or multipolygons.
        The faces within the polygons are removed.
    """
    logger.info('Clipping 2D mesh by polygon.')

    # Get xnodes, ynodes and edge_nodes
    xnodes = self.meshgeom.get_values('nodex', as_array=True)
    ynodes = self.meshgeom.get_values('nodey', as_array=True)
    edge_nodes = self.meshgeom.get_values('edge_nodes', as_array=True)

    # Clip
    logger.info('Clipping nodes.')
    xnodes, ynodes, edge_nodes = self.clip_nodes(xnodes, ynodes, edge_nodes, polygon, keep_outside=outside)

    # Update dimensions
    self.meshgeomdim.numnode = len(xnodes)
    self.meshgeomdim.numedge = len(edge_nodes)

    # Add nodes and links
    self.meshgeom.set_values('nodex', xnodes)
    self.meshgeom.set_values('nodey', ynodes)
    self.meshgeom.set_values('edge_nodes', np.ravel(edge_nodes).tolist())

    # Determine what the cells are
    logger.info('Finding cell faces within clipped nodes.')
    self._find_cells(self.meshgeom, dll_path=self.dll_path)

    # Clean nodes, this function deletes nodes based on no longer existing cells
    face_nodes = self.meshgeom.get_values('face_nodes', as_array=True)
    logger.info('Removing incomplete faces.')
    cleaned = self.clean_nodes(xnodes, ynodes, edge_nodes, face_nodes)

    # If cleaning leaves no whole cells, return None
    if cleaned is None:
        raise ValueError('Clipping the grid does not leave any nodes.')
    # Else, get the arguments from the returned tuple
    else:
        xnodes, ynodes, edge_nodes, face_nodes = cleaned

    logger.info('Update mesh with clipped nodes, edges and faces.')
    # Update dimensions
    self.meshgeomdim.numnode = len(xnodes)
    self.meshgeomdim.numedge = len(edge_nodes)

    # Add nodes and links
    self.meshgeom.set_values('nodex', xnodes)
    self.meshgeom.set_values('nodey', ynodes)
    self.meshgeom.set_values('edge_nodes', np.ravel(edge_nodes).tolist())
    self.meshgeom.set_values('face_nodes', np.ravel(face_nodes).tolist())

clip_nodes(xnodes, ynodes, edge_nodes, polygon, keep_outside=False)

Method that returns a set of nodes that is clipped within a given polygon.

Parameters

xnodes : X-coordinates of mesh nodes ynodes : Y-coordinates of mesh nodes edge_nodes : node ids of edges. These are assumed in 1-based numbering polygon : Polygon, Multipolygon or list of

Source code in hydrolib/dhydamo/geometry/mesh2d_gridgeom.py

def clip_nodes(self, xnodes, ynodes, edge_nodes, polygon, keep_outside=False):
    """
    Method that returns a set of nodes that is clipped within a given polygon.

    Parameters
    ----------
    xnodes : 
        X-coordinates of mesh nodes
    ynodes : 
        Y-coordinates of mesh nodes
    edge_nodes : 
        node ids of edges. These are assumed in 1-based numbering
    polygon : Polygon, Multipolygon or list of

    """

    # Find nodes within one of the polygons
    logger.info('Selecting nodes within polygon.')
    index = np.zeros(len(xnodes), dtype=bool)
    for poly in geometry.as_polygon_list(polygon):
        index = index | geometry.points_in_polygon(np.c_[xnodes, ynodes], poly)
    if keep_outside:
        index = ~index

    # Filter edge nodes, keep the ones that are within a polygon with both points
    nodeids_in_polygon = np.where(index)[0] + 1
    edge_nodes = edge_nodes[np.isin(edge_nodes, nodeids_in_polygon).all(axis=1)]

    # Remove edges that intersect the clipping polygons. This is needed for when the clipping polygon is narrow,
    # So it crosses an edge without containing one of the edge nodes
    logger.info('Creating LineString for each edge.')
    edges = [LineString(line) for line in np.c_[xnodes, ynodes][edge_nodes - 1]]

    logger.info('Checking for intersections.')
    isect = np.zeros(len(edge_nodes), dtype=bool)
    for poly in geometry.as_polygon_list(polygon):
        prepped = prep(poly.exterior)
        isect = isect | np.array([prepped.intersects(e) for e in edges])
        for interior in poly.interiors:
            prepped = prep(interior)
            isect = isect | np.array([prepped.intersects(e) for e in edges])
    edge_nodes = edge_nodes[~isect]

    # Remove edges that have a 1 count, until none are left
    logger.info('Remove edges with only a single node within the clip area.')
    unique, count = np.unique(edge_nodes, return_counts=True)
    while any(count == 1):
        count1nodes = unique[count == 1]
        edge_nodes = edge_nodes[~np.isin(edge_nodes, count1nodes).any(axis=1)]
        unique, count = np.unique(edge_nodes, return_counts=True)

    # Select remaining nodes
    xnodes = xnodes[unique - 1]
    ynodes = ynodes[unique - 1]

    # Make mapping for new id's
    new_id_mapping = {old_id: new_id + 1 for new_id, old_id in enumerate(unique)}
    edge_nodes = np.reshape([new_id_mapping[old_id] for old_id in edge_nodes.ravel()], edge_nodes.shape)

    return xnodes, ynodes, edge_nodes

faces_to_centroid()

The find_cells function does not always return face coordinates that are all accepted by the interactor. The coordinates are placed on the circumcenters, and placed on the face border in case of a coordinate that is outside the cells.

This function shifts the face coordinates towards the centroid (center of gravity) of the cell.

Source code in hydrolib/dhydamo/geometry/mesh2d_gridgeom.py

def faces_to_centroid(self):
    """
    The find_cells function does not always return face coordinates
    that are all accepted by the interactor. The coordinates
    are placed on the circumcenters, and placed on the face border
    in case of a coordinate that is outside the cells.

    This function shifts the face coordinates towards the centroid
    (center of gravity) of the cell.
    """
    # Calculate centroids
    centroids = np.vstack([cell.mean(axis=0) for cell in self.meshgeom.get_faces()]).T

    # Set values back to geometry
    self.meshgeom.set_values('facex', centroids[0])
    self.meshgeom.set_values('facey', centroids[1])

Rectangular

Bases: Mesh2D_GG

Source code in hydrolib/dhydamo/geometry/mesh2d_gridgeom.py

class Rectangular(Mesh2D_GG):

    def __init__(self, dll_path=None):
        Mesh2D_GG.__init__(self,dll_path=dll_path)
        self.meshgeomdim.maxnumfacenodes = 4
        self.dll_path=dll_path
        self.rotated = False

    def generate_grid(self, x0, y0, dx, dy, ncols, nrows, clipgeo=None, rotation=0):
        """
        Generate rectangular grid based on the origin (x0, y0) the cell sizes (dx, dy),
        the number of columns and rows (ncols, nrows) and a rotation in degrees (default=0)
        A geometry (clipgeo) can be given to clip the grid.
        """

        # Generate x and y spacing
        x = np.linspace(x0, x0 + dx * (ncols), ncols+1)
        y = np.linspace(y0, y0 + dy * (nrows), nrows+1)

        # Get nodes
        xnodes, ynodes = np.meshgrid(x, y)

        # Rotate
        if rotation != 0:
            self.rotated = True
            xnodes, ynodes = geometry.rotate_coordinates((x0, y0), np.radians(rotation), xnodes, ynodes)

        # Get nodes as list
        nodes = {crd: i+1 for i, crd in enumerate(zip(xnodes.ravel(), ynodes.ravel()))}

        # Create segments
        segments = list(zip(
            zip(*np.c_[xnodes[:, :-1].ravel(), ynodes[:, :-1].ravel()].T),
            zip(*np.c_[xnodes[:, 1:].ravel(), ynodes[:, 1:].ravel()].T)
        ))
        segments += list(zip(
            zip(*np.c_[xnodes[1:, :].ravel(), ynodes[1:, :].ravel()].T),
            zip(*np.c_[xnodes[:-1, :].ravel(), ynodes[:-1, :].ravel()].T)
        ))

        # Get links
        edge_nodes = np.asarray([(nodes[s[0]], nodes[s[1]]) for s in segments])

        # Rvael the nodes
        xnodes = xnodes.ravel()
        ynodes = ynodes.ravel()

        # Add to mesh
        dimensions = meshgeomdim()
        dimensions.dim = 2
        geometries = meshgeom(dimensions)

        # Clip
        if clipgeo is not None:
            xnodes, ynodes, edge_nodes = self.clip_nodes(xnodes, ynodes, edge_nodes, clipgeo)

        # Update dimensions
        dimensions.numnode = len(xnodes)
        dimensions.numedge = len(edge_nodes)

        # Add nodes and links
        geometries.set_values('nodex', xnodes)
        geometries.set_values('nodey', ynodes)
        geometries.set_values('edge_nodes', np.ravel(edge_nodes).tolist())

        # Determine what the cells are
        self._find_cells(geometries, dll_path=self.dll_path)

        # Clip
        if clipgeo is not None:
            # Clean nodes, this function deletes nodes based on no longer existing cells
            face_nodes = geometries.get_values('face_nodes', as_array=True)
            cleaned = self.clean_nodes(xnodes, ynodes, edge_nodes, face_nodes)
            # If cleaning leaves no whole cells, return None
            if cleaned is None:
                return None
            # Else, get the arguments from the returned tuple
            else:
                xnodes, ynodes, edge_nodes, face_nodes = cleaned

        # Update dimensions
        dimensions.numnode = len(xnodes)
        dimensions.numedge = len(edge_nodes)
        dimensions.maxnumfacenodes = 4

        # Add nodes and links
        geometries.set_values('nodex', xnodes)
        geometries.set_values('nodey', ynodes)
        geometries.set_values('edge_nodes', np.ravel(edge_nodes).tolist())
        geometries.set_values('face_nodes', np.ravel(face_nodes).tolist())

        # Add to mesh
        self.meshgeom.add_from_other(geometries)

    def generate_within_polygon(self, polygon, cellsize, rotation=0):
        """
        Function to generate a grid within a polygon. It uses the function
        'generate_grid' but automatically detects the extent.

        Parameters
        ----------
        polygon : (list of) shapely.geometry.Polygon or a shapely.geometry.MultiPolygon
        """

        polygons = geometry.as_polygon_list(polygon)

        logger.info(
            "Generating grid with cellsize %s m and rotation %s degrees.",
            cellsize,
            rotation,
        )
        convex = MultiPolygon(polygons).convex_hull

        # In case of a rotation, extend the grid far enough to make sure
        if rotation != 0:
            # Find a box that contains the whole polygon
            (x0, y0), xsize, ysize = geometry.minimum_bounds_fixed_rotation(convex, rotation)
        else:
            xsize = convex.bounds[2] - convex.bounds[0]
            ysize= convex.bounds[3] - convex.bounds[1]
            x0, y0 = convex.bounds[0], convex.bounds[1]

        self.generate_grid(
            x0=x0,
            y0=y0,
            dx=cellsize,
            dy=cellsize,
            ncols=int(xsize / cellsize) + 1,
            nrows=int(ysize / cellsize) + 1,
            clipgeo=polygon,
            rotation=rotation
        )

    def refine(self, polygon, level, cellsize, debug=False, dflowfm_path=None):
        """
        Parameters
        ----------
        polygon : list of tuples
            Polygon in which to refine
        level : int
            Number of times to refine
        """

        if self.rotated:
            raise NotImplementedError('Mesh refinement does not work for rotation grids.')

        # checks.check_argument(polygon, 'polygon', (list, Polygon, MultiPolygon))

        # Save network as grid
        temppath = 'tempgrid_net.nc'
        gridio.to_netcdf_old(self.meshgeom, temppath)

        if isinstance(level, int):
            level = [level]
            polygon = [polygon]

        factor = 2 ** max(level)
        stepsize = cellsize / factor
        dxmin = stepsize
        dtmax = stepsize / (9.81) ** 0.5


        xnodes = self.meshgeom.get_values('nodex')
        ynodes = self.meshgeom.get_values('nodey')
        x0 = min(xnodes)
        ncols = int((max(xnodes) - x0) / cellsize)
        y0 = min(ynodes)
        nrows = int((max(ynodes) - y0) / cellsize)

        arr = np.ones((nrows * factor, ncols * factor)) * 4**max(level)

        for poly, lvl in zip(polygon, level):
            mask = geometry.geometry_to_mask(poly, (x0, y0), stepsize, shape=(nrows * factor, ncols * factor))
            arr[mask] = 4 ** (max(level) - lvl)

        # Create ascii grid based on polygon
        header = (
            'nCols        {0:d}\n'
            'nRows        {1:d}\n'
            'xllCorner    {2:.6f}\n'
            'yllCorner    {3:.6f}\n'
            'CellSize     {4:.6f}\n'
            'nodata_value 0\n'
        ).format(ncols * factor, nrows * factor, x0, y0 + stepsize, stepsize)

        values = '\n'.join(' '.join('%d' %x for x in y) for y in arr.squeeze())

        with open('temp_ascgrid.asc', 'w') as f:
            f.write(header + values)

        # Get dflow fm on system path
        if dflowfm_path is None or not os.path.exists(dflowfm_path):
            raise OSError('dflowfm.exe not provided or not found.')        

        # Construct statement        
        statement = f'"{dflowfm_path}" --refine:hmin={dxmin}:dtmax={dtmax}:connect=1:outsidecell=1 {temppath} temp_ascgrid.asc'
        if debug:
            statement += ' > log.txt'
        # Execute statement
        os.system(statement)

        if not debug:
            os.remove('temp_ascgrid.asc')

        # Read geometry and dimensions again
        if not os.path.exists('out_net.nc'):
            raise OSError(
                (f'Could not find "out_net.nc" in directory "{os.getcwd()}". '
                 f'The file should have been generated with the statement: "{statement}". '
                 'You can test this manually by running this statement from a command prompt launched from the directory named in this message.'))
        gridio.from_netcdf_old(self.meshgeom, 'out_net.nc')

        # Find cells
        self._find_cells(self.meshgeom, dll_path=self.dll_path)

        os.remove('out_net.nc')
        os.remove(temppath)

generate_grid(x0, y0, dx, dy, ncols, nrows, clipgeo=None, rotation=0)

Generate rectangular grid based on the origin (x0, y0) the cell sizes (dx, dy), the number of columns and rows (ncols, nrows) and a rotation in degrees (default=0) A geometry (clipgeo) can be given to clip the grid.

Source code in hydrolib/dhydamo/geometry/mesh2d_gridgeom.py

def generate_grid(self, x0, y0, dx, dy, ncols, nrows, clipgeo=None, rotation=0):
    """
    Generate rectangular grid based on the origin (x0, y0) the cell sizes (dx, dy),
    the number of columns and rows (ncols, nrows) and a rotation in degrees (default=0)
    A geometry (clipgeo) can be given to clip the grid.
    """

    # Generate x and y spacing
    x = np.linspace(x0, x0 + dx * (ncols), ncols+1)
    y = np.linspace(y0, y0 + dy * (nrows), nrows+1)

    # Get nodes
    xnodes, ynodes = np.meshgrid(x, y)

    # Rotate
    if rotation != 0:
        self.rotated = True
        xnodes, ynodes = geometry.rotate_coordinates((x0, y0), np.radians(rotation), xnodes, ynodes)

    # Get nodes as list
    nodes = {crd: i+1 for i, crd in enumerate(zip(xnodes.ravel(), ynodes.ravel()))}

    # Create segments
    segments = list(zip(
        zip(*np.c_[xnodes[:, :-1].ravel(), ynodes[:, :-1].ravel()].T),
        zip(*np.c_[xnodes[:, 1:].ravel(), ynodes[:, 1:].ravel()].T)
    ))
    segments += list(zip(
        zip(*np.c_[xnodes[1:, :].ravel(), ynodes[1:, :].ravel()].T),
        zip(*np.c_[xnodes[:-1, :].ravel(), ynodes[:-1, :].ravel()].T)
    ))

    # Get links
    edge_nodes = np.asarray([(nodes[s[0]], nodes[s[1]]) for s in segments])

    # Rvael the nodes
    xnodes = xnodes.ravel()
    ynodes = ynodes.ravel()

    # Add to mesh
    dimensions = meshgeomdim()
    dimensions.dim = 2
    geometries = meshgeom(dimensions)

    # Clip
    if clipgeo is not None:
        xnodes, ynodes, edge_nodes = self.clip_nodes(xnodes, ynodes, edge_nodes, clipgeo)

    # Update dimensions
    dimensions.numnode = len(xnodes)
    dimensions.numedge = len(edge_nodes)

    # Add nodes and links
    geometries.set_values('nodex', xnodes)
    geometries.set_values('nodey', ynodes)
    geometries.set_values('edge_nodes', np.ravel(edge_nodes).tolist())

    # Determine what the cells are
    self._find_cells(geometries, dll_path=self.dll_path)

    # Clip
    if clipgeo is not None:
        # Clean nodes, this function deletes nodes based on no longer existing cells
        face_nodes = geometries.get_values('face_nodes', as_array=True)
        cleaned = self.clean_nodes(xnodes, ynodes, edge_nodes, face_nodes)
        # If cleaning leaves no whole cells, return None
        if cleaned is None:
            return None
        # Else, get the arguments from the returned tuple
        else:
            xnodes, ynodes, edge_nodes, face_nodes = cleaned

    # Update dimensions
    dimensions.numnode = len(xnodes)
    dimensions.numedge = len(edge_nodes)
    dimensions.maxnumfacenodes = 4

    # Add nodes and links
    geometries.set_values('nodex', xnodes)
    geometries.set_values('nodey', ynodes)
    geometries.set_values('edge_nodes', np.ravel(edge_nodes).tolist())
    geometries.set_values('face_nodes', np.ravel(face_nodes).tolist())

    # Add to mesh
    self.meshgeom.add_from_other(geometries)

generate_within_polygon(polygon, cellsize, rotation=0)

Function to generate a grid within a polygon. It uses the function 'generate_grid' but automatically detects the extent.

Parameters

polygon : (list of) shapely.geometry.Polygon or a shapely.geometry.MultiPolygon

Source code in hydrolib/dhydamo/geometry/mesh2d_gridgeom.py

def generate_within_polygon(self, polygon, cellsize, rotation=0):
    """
    Function to generate a grid within a polygon. It uses the function
    'generate_grid' but automatically detects the extent.

    Parameters
    ----------
    polygon : (list of) shapely.geometry.Polygon or a shapely.geometry.MultiPolygon
    """

    polygons = geometry.as_polygon_list(polygon)

    logger.info(
        "Generating grid with cellsize %s m and rotation %s degrees.",
        cellsize,
        rotation,
    )
    convex = MultiPolygon(polygons).convex_hull

    # In case of a rotation, extend the grid far enough to make sure
    if rotation != 0:
        # Find a box that contains the whole polygon
        (x0, y0), xsize, ysize = geometry.minimum_bounds_fixed_rotation(convex, rotation)
    else:
        xsize = convex.bounds[2] - convex.bounds[0]
        ysize= convex.bounds[3] - convex.bounds[1]
        x0, y0 = convex.bounds[0], convex.bounds[1]

    self.generate_grid(
        x0=x0,
        y0=y0,
        dx=cellsize,
        dy=cellsize,
        ncols=int(xsize / cellsize) + 1,
        nrows=int(ysize / cellsize) + 1,
        clipgeo=polygon,
        rotation=rotation
    )

refine(polygon, level, cellsize, debug=False, dflowfm_path=None)

Parameters

polygon : list of tuples Polygon in which to refine level : int Number of times to refine

Source code in hydrolib/dhydamo/geometry/mesh2d_gridgeom.py

def refine(self, polygon, level, cellsize, debug=False, dflowfm_path=None):
    """
    Parameters
    ----------
    polygon : list of tuples
        Polygon in which to refine
    level : int
        Number of times to refine
    """

    if self.rotated:
        raise NotImplementedError('Mesh refinement does not work for rotation grids.')

    # checks.check_argument(polygon, 'polygon', (list, Polygon, MultiPolygon))

    # Save network as grid
    temppath = 'tempgrid_net.nc'
    gridio.to_netcdf_old(self.meshgeom, temppath)

    if isinstance(level, int):
        level = [level]
        polygon = [polygon]

    factor = 2 ** max(level)
    stepsize = cellsize / factor
    dxmin = stepsize
    dtmax = stepsize / (9.81) ** 0.5


    xnodes = self.meshgeom.get_values('nodex')
    ynodes = self.meshgeom.get_values('nodey')
    x0 = min(xnodes)
    ncols = int((max(xnodes) - x0) / cellsize)
    y0 = min(ynodes)
    nrows = int((max(ynodes) - y0) / cellsize)

    arr = np.ones((nrows * factor, ncols * factor)) * 4**max(level)

    for poly, lvl in zip(polygon, level):
        mask = geometry.geometry_to_mask(poly, (x0, y0), stepsize, shape=(nrows * factor, ncols * factor))
        arr[mask] = 4 ** (max(level) - lvl)

    # Create ascii grid based on polygon
    header = (
        'nCols        {0:d}\n'
        'nRows        {1:d}\n'
        'xllCorner    {2:.6f}\n'
        'yllCorner    {3:.6f}\n'
        'CellSize     {4:.6f}\n'
        'nodata_value 0\n'
    ).format(ncols * factor, nrows * factor, x0, y0 + stepsize, stepsize)

    values = '\n'.join(' '.join('%d' %x for x in y) for y in arr.squeeze())

    with open('temp_ascgrid.asc', 'w') as f:
        f.write(header + values)

    # Get dflow fm on system path
    if dflowfm_path is None or not os.path.exists(dflowfm_path):
        raise OSError('dflowfm.exe not provided or not found.')        

    # Construct statement        
    statement = f'"{dflowfm_path}" --refine:hmin={dxmin}:dtmax={dtmax}:connect=1:outsidecell=1 {temppath} temp_ascgrid.asc'
    if debug:
        statement += ' > log.txt'
    # Execute statement
    os.system(statement)

    if not debug:
        os.remove('temp_ascgrid.asc')

    # Read geometry and dimensions again
    if not os.path.exists('out_net.nc'):
        raise OSError(
            (f'Could not find "out_net.nc" in directory "{os.getcwd()}". '
             f'The file should have been generated with the statement: "{statement}". '
             'You can test this manually by running this statement from a command prompt launched from the directory named in this message.'))
    gridio.from_netcdf_old(self.meshgeom, 'out_net.nc')

    # Find cells
    self._find_cells(self.meshgeom, dll_path=self.dll_path)

    os.remove('out_net.nc')
    os.remove(temppath)

models

GeometryList

Bases: GeometryList

Source code in hydrolib/dhydamo/geometry/models.py

class GeometryList(GeometryListMK):
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

    @classmethod
    def from_geometry(cls, geometry):
        if isinstance(geometry, Polygon):
            return cls.from_polygon(geometry)
        elif isinstance(geometry, MultiPolygon):
            return cls.from_multipolygon(geometry)
        elif isinstance(geometry, LineString):
            return cls.from_linestring(geometry)
        elif isinstance(geometry, MultiLineString):
            return cls.from_multilinestring(geometry)
        elif isinstance(geometry, Point):
            return cls.from_point(geometry)
        elif isinstance(geometry, MultiPoint):
            return cls.from_multipoint(geometry)
        else:
            raise TypeError(f"Geometry type {type(geometry)} not understood.")

    @classmethod
    def _from_simple(cls, geometry: LineString | Point):
        # Extract coordinates from geometry
        x_crds, y_crds = np.array(geometry.coords[:]).T
        gl = cls(x_coordinates=x_crds, y_coordinates=y_crds)
        return gl

    @classmethod
    def from_linestring(cls, linestring: LineString):
        return cls._from_simple(linestring)

    @classmethod
    def from_point(cls, point: Point):
        return cls._from_simple(point)

    @classmethod
    def from_polygon(cls, polygon: Polygon):
        # Create a list of coordinate lists
        # Add exterior
        x_ext, y_ext = np.array(polygon.exterior.coords[:]).T
        x_crds = [x_ext]
        y_crds = [y_ext]


        if not hasattr(cls, 'inner_outer_separator'):
            cls.inner_outer_separator = -998

        # Add interiors, seperated by inner_outer_separator
        for interior in polygon.interiors:
            x_int, y_int = np.array(interior.coords[:]).T
            x_crds.append([cls.inner_outer_separator])
            x_crds.append(x_int)
            y_crds.append([cls.inner_outer_separator])
            y_crds.append(y_int)
        gl = cls(
            x_coordinates=np.concatenate(x_crds), y_coordinates=np.concatenate(y_crds)
        )
        return gl

    @classmethod
    def _from_multigeometry(cls, multigeometry: MultiLineString | MultiPolygon):
        x_crds = []
        y_crds = []
        cls.inner_outer_separator = -998
        cls.geometry_separator = -999
        # Create a GeometryList for every polygon in the multipolygon
        for i, geometry in enumerate(multigeometry.geoms):
            gl = cls.from_geometry(geometry)
            # Add geometry seperator for every geometrylist, except the first one
            if i != 0:
                x_crds.append([cls.geometry_separator])
                y_crds.append([cls.geometry_separator])
            # Add coordinates of polygon GeometryList
            x_crds.append(gl.x_coordinates)
            y_crds.append(gl.y_coordinates)

        gl = cls(
            x_coordinates=np.concatenate(x_crds), y_coordinates=np.concatenate(y_crds)
        )
        return gl

    @classmethod
    def from_multipolygon(cls, multipolygon: MultiPolygon):
        return cls._from_multigeometry(multipolygon)

    @classmethod
    def from_multilinestring(cls, multilinestring: LineString):
        return cls._from_multigeometry(multilinestring)

    @classmethod
    def from_multipoint(cls, multipoint: Point):
        return cls._from_multigeometry(multipoint)

    def _to_polygon(
        self, geometries: list[GeometryListMK], is_multi: bool
    ) -> Polygon | MultiPolygon:
        polygons = []
        for geometry in geometries:
            parts = [
                np.stack([p.x_coordinates, p.y_coordinates], axis=1)
                for p in split_by(geometry, self.inner_outer_separator)
            ]
            polygons.append(Polygon(shell=parts[0], holes=parts[1:]))
        if is_multi:
            return MultiPolygon(polygons)
        else:
            return polygons[0] if len(polygons)==1 else None

    def _to_linestring(
        self, geometries: list[GeometryListMK], is_multi: bool
    ) -> LineString | MultiLineString:
        linestrings = [
            LineString(np.stack([p.x_coordinates, p.y_coordinates], axis=1))
            for p in geometries
        ]
        if is_multi:
            return MultiLineString(linestrings)
        else:
            return linestrings[0]

    def _to_points(
        self, geometries: list[GeometryListMK], is_multi: bool
    ) -> Point | MultiPoint:
        points = [
            Point(np.stack([p.x_coordinates, p.y_coordinates], axis=1))
            for p in geometries
        ]
        if is_multi:
            return MultiPoint(points)
        else:
            return points[0]

    def to_geometry(self):
        geometries = [
            geo
            for geo in split_by(self, self.geometry_separator)
            if geo.x_coordinates.size > 0
        ]
        is_multi = len(geometries) > 1

        for geometry_list in geometries:
            # Check if polygon, by comparing first and last coordinates
            exterior = split_by(geometry_list, self.inner_outer_separator)[0]
            is_polygon = exterior.x_coordinates[0] == exterior.x_coordinates[-1]
            # Check if linestring, by checking the length of coordinate sequences
            is_linestring = (len(geometry_list.x_coordinates) > 1) and (not is_polygon)

            if is_polygon:
                return self._to_polygon(geometries, is_multi)

            elif is_linestring:
                return self._to_linestring(geometries, is_multi)

            else:
                return self._to_point(geometries, is_multi)

    @property
    def geoms(self):
        """Returns GeometryList objects based on spliited by geometries"""
        for geometrylist in split_by(self, self.geometry_separator):
            # yield as shapely geometry
            yield GeometryList(**geometrylist.__dict__).to_geometry()

geoms property

Returns GeometryList objects based on spliited by geometries

rasterstats

check_geodateframe_rasterstats(facedata)

Check for type, columns and coordinates

Source code in hydrolib/dhydamo/geometry/rasterstats.py

def check_geodateframe_rasterstats(facedata):
    """
    Check for type, columns and coordinates
    """
    if not isinstance(facedata, gpd.GeoDataFrame):
        raise TypeError("facedata should be type GeoDataFrame")

    # Check if facedata has required columns
    if ("facex" not in facedata.columns) or ("facey" not in facedata.columns):
        xy = list(zip(*[pt.coords[0] for pt in facedata.geometry.centroid]))
        facedata["facex"] = xy[0]
        facedata["facey"] = xy[1]

    # Check if coordinates are present.
    if "crds" not in facedata.columns:
        facedata["crds"] = [row.coords[:] for row in facedata.geometry]

compress(path)

Function re-save an existing raster file with compression.

Parameters

path : str Path to raster file. File is overwritten with compress variant.

Source code in hydrolib/dhydamo/geometry/rasterstats.py

def compress(path):
    """
    Function re-save an existing raster file with compression.

    Parameters
    ----------
    path : str
        Path to raster file. File is overwritten with compress variant.
    """
    # Compress
    with rasterio.open(path, "r") as f:
        arr = f.read()
        out_meta = f.meta.copy()
        out_meta["compress"] = "deflate"
    with rasterio.open(path, "w", **out_meta) as f:
        f.write(arr)

raster_in_parts(f: rasterio.io.DatasetReader, ncols: int, nrows: int, facedata) -> RasterPart

Certain rasters are too big to read into memory at once. This function helps splitting them in equal parts of (+- ncols x nrows pixels)

If facedata is given, each part is extended such that whole faces are covered by the parts

Parameters:

Name	Type	Description	Default
f	_type_	description	required
ncols	_type_	description	required
nrows	_type_	description	required
facedata	_type_	description	required

Yields:

Name	Type	Description
_type_	RasterPart	description

Source code in hydrolib/dhydamo/geometry/rasterstats.py

def raster_in_parts(
    f: rasterio.io.DatasetReader, ncols: int, nrows: int, facedata
) -> RasterPart:
    """Certain rasters are too big to read into memory at once.
    This function helps splitting them in equal parts of (+- ncols x nrows pixels)

    If facedata is given, each part is extended such that whole faces
    are covered by the parts

    Args:
        f (_type_): _description_
        ncols (_type_): _description_
        nrows (_type_): _description_
        facedata (_type_): _description_

    Yields:
        _type_: _description_
    """
    nx = max(1, f.shape[1] // ncols)
    ny = max(1, f.shape[0] // nrows)

    xparts = np.linspace(0, f.shape[1], nx + 1).astype(int)
    yparts = np.linspace(0, f.shape[0], ny + 1).astype(int)

    pts = facedata[["facex", "facey"]].to_numpy()

    for ix, iy in product(range(nx), range(ny)):
        part = RasterPart(
            f,
            xmin=xparts[ix],
            ymin=yparts[iy],
            xmax=xparts[ix + 1],
            ymax=yparts[iy + 1],
        )

        if facedata is not None:
            # For each part, get points in part.
            # For narrow/diagonal shapes the points in part can be limited
            idx = part.get_pts_in_part(pts)
            if not idx.any():
                continue

            crds = facedata["crds"].tolist()
            ll = list(zip(*[crds[i].min(axis=0) for i in np.where(idx)[0]]))
            ur = list(zip(*[crds[i].max(axis=0) for i in np.where(idx)[0]]))
            bounds = (min(ll[0]), min(ll[1]), max(ur[0]), max(ur[1]))

            # Get new part based on extended bounds
            part = RasterPart.from_bounds(f, bounds)

            # Add the cell centers within the window as index to the part
            part.idx = idx

        yield part

raster_stats_fine_cells(rasterpath: str | Path, facedata, stats=['mean'])

Calculate statistic from a raster, where the raster resoltion is (much) smaller than the cell size.

Parameters

rasterpath : str Path to raster file facedata : geopandas.GeoDataFrame Dataframe with polygons in which the raster statistics are derived. stats : list List of statistics to retrieve. Should be numpy functions that require one argument

Source code in hydrolib/dhydamo/geometry/rasterstats.py

def raster_stats_fine_cells(rasterpath: str | Path, facedata, stats=["mean"]):
    """
    Calculate statistic from a raster, where the raster resoltion is (much)
    smaller than the cell size.

    Parameters
    ----------
    rasterpath : str
        Path to raster file
    facedata : geopandas.GeoDataFrame
        Dataframe with polygons in which the raster statistics are derived.
    stats : list
        List of statistics to retrieve. Should be numpy functions that require one argument
    """

    # Create empty array for stats
    stat_array = {stat: {} for stat in stats + ["count"]}

    # Check geometries
    check_geodateframe_rasterstats(facedata)

    # Open raster file
    with rasterio.open(rasterpath, "r") as f:
        # Split file in parts based on shape
        parts = raster_in_parts(f, ncols=250, nrows=250, facedata=facedata)

        for prt in parts:
            # Get values from array
            arr = prt.read(1)
            valid = arr != f.nodata
            if not valid.any():
                continue

            # Rasterize the cells in the part
            cellidx_sel = rasterize_cells(facedata.loc[prt.idx], prt)
            assert cellidx_sel.shape == valid.shape
            cellidx_sel[~valid] = 0
            valid = cellidx_sel != 0

            for cell_idx in np.unique(cellidx_sel):
                if cell_idx == 0:
                    continue

                # Mask
                cellmask = cellidx_sel == cell_idx
                if not cellmask.any():
                    continue

                # Get bottom values
                bottom = arr[cellmask]

                # For each statistic, get the values and add
                for stat in stats:
                    stat_array[stat][cell_idx] = getattr(np, stat)(bottom)

                stat_array["count"][cell_idx] = len(bottom)

        # Cast to pandas dataframe
        df = pd.DataFrame.from_dict(stat_array).reindex(index=facedata.index)

    return df

waterdepth_ahn(dempath, facedata, outpath, column)

Function that combines a dem and water levels to a water depth raster. No sub grid correction is done.

Parameters

dempath : str Path to raster file with terrain level facedata : gpd.GeoDataFrame GeoDataFrame with at least the cell geometry and a column with water levels outpath : str Path to output raster file column : str Name of the column with the water level data

Source code in hydrolib/dhydamo/geometry/rasterstats.py

def waterdepth_ahn(dempath, facedata, outpath, column):
    """
    Function that combines a dem and water levels to a water
    depth raster. No sub grid correction is done.

    Parameters
    ----------
    dempath : str
        Path to raster file with terrain level
    facedata : gpd.GeoDataFrame
        GeoDataFrame with at least the cell geometry and
        a column with water levels
    outpath : str
        Path to output raster file
    column : str
        Name of the column with the water level data
    """

    # Open raster file
    with rasterio.open(dempath, "r") as f:
        first = True
        out_meta = f.meta.copy()

        # Split file in parts based on shape
        parts = raster_in_parts(f, ncols=250, nrows=250, facedata=facedata)

        for prt in parts:
            # Get values from array
            arr = prt.read(1)
            valid = arr != f.nodata
            if not valid.any():
                continue

            cellidx_sel = rasterize_cells(facedata.loc[prt.idx], prt)
            cellidx_sel[~valid] = 0
            valid = cellidx_sel != 0

            # Create array to assign water levels
            wlev_subgr = np.zeros(cellidx_sel.shape, dtype=out_meta["dtype"])

            for cell_idx in np.unique(cellidx_sel):
                if cell_idx == 0:
                    continue
                # Mask
                cellmask = cellidx_sel == cell_idx
                if not cellmask.any():
                    continue
                # Add values for cell to raster
                wlev_subgr[cellmask] = facedata.at[cell_idx, column]

            # Write to output raster
            with rasterio.open(outpath, "w" if first else "r+", **out_meta) as dst:
                # Determine water depth
                if not first:
                    wdep_subgr = dst.read(1, window=prt.window)
                else:
                    wdep_subgr = np.ones_like(wlev_subgr) * f.nodata
                    first = False
                wdep_subgr[valid] = np.maximum(wlev_subgr - arr, 0).astype(
                    out_meta["dtype"]
                )[valid]
                dst.write(wdep_subgr[None, :, :], window=prt.window)

    compress(outpath)

spatial

find_nearest_branch(branches, geometries, method='overal', maxdist=5)

Method to determine nearest branch for each geometry. The nearest branch can be found by finding t from both ends (ends) or the nearest branch from the geometry as a whole (overal), the centroid (centroid), or intersecting (intersect).

Parameters

branches : geopandas.GeoDataFrame Geodataframe with branches geometries : geopandas.GeoDataFrame Geodataframe with geometries to snap method='overal' : str Method for determine branch maxdist=5 : int or float Maximum distance for finding nearest geometry minoffset : int or float Minimum offset from the end of the corresponding branch in case of method=equal

Source code in hydrolib/dhydamo/geometry/spatial.py

def find_nearest_branch(branches, geometries, method='overal', maxdist=5):
    """
    Method to determine nearest branch for each geometry.
    The nearest branch can be found by finding t from both ends (ends) or the nearest branch from the geometry
    as a whole (overal), the centroid (centroid), or intersecting (intersect).

    Parameters
    ----------
    branches : geopandas.GeoDataFrame
        Geodataframe with branches
    geometries : geopandas.GeoDataFrame
        Geodataframe with geometries to snap
    method='overal' : str
        Method for determine branch
    maxdist=5 : int or float
        Maximum distance for finding nearest geometry
    minoffset : int or float
        Minimum offset from the end of the corresponding branch in case of method=equal
    """
    # Check if method is in allowed methods
    allowed_methods = ['intersecting', 'overal', 'centroid', 'ends']
    if method not in allowed_methods:
        raise NotImplementedError(f'Method "{method}" not implemented.')

    # Add columns if not present
    if 'branch_id' not in geometries.columns:
        geometries['branch_id'] = ''
    if 'branch_offset' not in geometries.columns:
        geometries['branch_offset'] = np.nan

    if method == 'intersecting':
        # Determine intersection geometries per branch
        geobounds = geometries.bounds.values.T
        for branch in branches.itertuples():
            selectie = geometries.loc[possibly_intersecting(geobounds, branch.geometry)].copy()
            intersecting = selectie.loc[selectie.intersects(branch.geometry).values]

            # For each geometrie, determine offset along branch
            for geometry in intersecting.itertuples():
                # Determine distance of profile line along branch
                geometries.at[geometry.Index, 'branch_id'] = branch.Index

                # Calculate offset
                branchgeo = branch.geometry
                mindist = min(0.1, branchgeo.length / 2.)
                offset = round(branchgeo.project(branchgeo.intersection(geometry.geometry).centroid), 3)
                offset = max(mindist, min(branchgeo.length - mindist, offset))
                geometries.at[geometry.Index, 'branch_offset'] = offset

    else:
        branch_bounds = branches.bounds.values.T
        # In case of looking for the nearest, it is easier to iteratie over the geometries instead of the branches
        for geometry in geometries.itertuples():
            # Find near branches
            nearidx = possibly_intersecting(branch_bounds, geometry.geometry, buffer=maxdist)
            selectie = branches.loc[nearidx]

            if method == 'overal':
                # Determine distances to branches
                dist = selectie.distance(geometry.geometry)
            elif method == 'centroid':
                # Determine distances to branches
                dist = selectie.distance(geometry.geometry.centroid)
            elif method == 'ends':
                # Since a culvert can cross a channel, it is
                crds = geometry.geometry.coords[:]
                dist = (selectie.distance(Point(*crds[0])) + selectie.distance(Point(*crds[-1]))) * 0.5

            # Determine nearest
            if dist.min() < maxdist:
                branchidxmin = dist.idxmin()
                geometries.at[geometry.Index, 'branch_id'] = branchidxmin
                if isinstance(geometry.geometry, Point):
                    geo = geometry.geometry
                else:
                    geo = geometry.geometry.centroid

                if (dist < maxdist).sum() > 1:
                    logger.info(
                        "Geometry centroid %s has multiple branches with maxdist=%s and method=%s: %s",
                        geo,
                        maxdist,
                        method,
                        dist.index[dist < maxdist].to_list(),
                    )

                # Calculate offset
                branchgeo = branches.at[branchidxmin, 'geometry']
                mindist = min(0.1, branchgeo.length / 2.)
                offset = max(mindist, min(branchgeo.length - mindist, round(branchgeo.project(geo), 3)))
                geometries.at[geometry.Index, 'branch_offset'] = offset

get_voronoi_around_nodes(nodes: np.ndarray, facedata: gpd.GeoDataFrame) -> gpd.GeoDataFrame

Creates voronoi polygons around face nodes.

Parameters:

Name	Type	Description	Default
nodes	ndarray	xy coordinates of face nodes	required
facedata	GeoDataFrame	GeoDataFrame with face properties	required

Returns:

Type	Description
GeoDataFrame	gpd.GeoDataFrame: Creating GeoDataFrame with created polygons and their properties

Source code in hydrolib/dhydamo/geometry/spatial.py

def get_voronoi_around_nodes(nodes: np.ndarray, facedata: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Creates voronoi polygons around face nodes.

    Args:
        nodes (np.ndarray): xy coordinates of face nodes
        facedata (gpd.GeoDataFrame): GeoDataFrame with face properties

    Returns:
        gpd.GeoDataFrame: Creating GeoDataFrame with created polygons and their properties
    """
    # Creat voronoi polygon
    # Add border to limit polygons
    border = box(nodes[:, 0].min(), nodes[:, 1].min(), nodes[:, 0].max(), nodes[:, 1].max()).buffer(1000).exterior
    borderpts = [border.interpolate(dist).coords[0] for dist in np.linspace(0, border.length, max(20, int(border.length / 100)))]
    vor = Voronoi(points=nodes.tolist()+borderpts)
    clippoly = facedata.union_all()
    # Get lines
    lines = []
    for poly in common.as_polygon_list(clippoly):
        lines.append(poly.exterior)
        lines.extend([line for line in poly.interiors])
    linesprep = prep(MultiLineString(lines))
    clipprep = prep(clippoly)

    # Collect polygons
    data = []
    for (pr, pt) in zip(vor.point_region, nodes):
        region = vor.regions[pr]
        if pr == -1:
            break
        while -1 in region:
            region.remove(-1)
        if len(region) < 3:
            continue
        crds = vor.vertices[region]
        if clipprep.intersects(Point(pt)):
            poly = Polygon(crds)
            if linesprep.intersects(poly):
                poly = poly.intersection(clippoly)
                if isinstance(poly, MultiPolygon):
                    poly = poly.buffer(0.001)
                if isinstance(poly, MultiPolygon):
                    logger.warning('Got multipolygon when clipping voronoi polygon. Only adding coordinates for largest of the polygons.')
                    poly = poly[np.argmax([p.area for p in common.as_polygon_list(poly)])]
                crds = np.vstack(poly.exterior.coords[:])
            data.append({'geometry': poly, 'crds': crds})

    # Limit to model extend
    facedata = gpd.GeoDataFrame(data)
    facedata.index=np.arange(len(nodes), dtype=np.uint32) + 1

    return facedata

minimum_bounds_fixed_rotation(polygon, angle)

Get the minimum box for a polygon with a given axes rotation.

Parameters

polygon : shapely.geometry.Polygon Polygon that is rotated angle : int or float Rotation of the polygon in degrees

Returns

tuple Tuple with origin (x, y), xsize and ysize

Source code in hydrolib/dhydamo/geometry/spatial.py

def minimum_bounds_fixed_rotation(polygon, angle):
    """Get the minimum box for a polygon with a given axes rotation.

    Parameters
    ----------
    polygon : shapely.geometry.Polygon
        Polygon that is rotated
    angle : int or float
        Rotation of the polygon in degrees

    Returns
    -------
    tuple
        Tuple with origin (x, y), xsize and ysize
    """
    # Determine spinning point
    spinpt = (polygon.envelope.bounds[0], polygon.envelope.bounds[1])

    # Rotate clip polygon with rotation, get envelope and rotate back.
    rotbox1 = affinity.rotate(polygon, angle=angle, origin=spinpt).envelope

    # Determine size of grid
    xsize = rotbox1.bounds[2] - rotbox1.bounds[0]
    ysize = rotbox1.bounds[3] - rotbox1.bounds[1]

    # Rotate again, and get origin
    rotbox2 = affinity.rotate(rotbox1, angle=-angle, origin=spinpt)
    origin = rotbox2.exterior.coords[0]

    return origin, xsize, ysize

orthogonal_line(line: LineString, offset: float, width: float = 1.0) -> list[tuple[float]]

Parameters

line : shapely.geometry.LineString Line geometry object on which the orthogonal line is drawn offset : float Offset of the orthogonal line along line width : float Width of the orthogonal line

Returns

line : list List with coordinate tuples

Source code in hydrolib/dhydamo/geometry/spatial.py

def orthogonal_line(line: LineString, offset: float, width: float=1.0) -> list[tuple[float]]:
    """
    Parameters
    ----------
    line : shapely.geometry.LineString
        Line geometry object on which the orthogonal line is drawn
    offset : float
        Offset of the orthogonal line along line
    width : float
        Width of the orthogonal line

    Returns
    -------
    line : list
        List with coordinate tuples
    """

    # Determine angle at offset
    angle = np.angle(complex(*np.diff([
        line.interpolate(offset - 0.001).coords[0][:2],
        line.interpolate(offset + 0.001).coords[0][:2]
    ], axis=0)[0])) + 0.5 * np.pi

    # Create new line
    pt = line.interpolate(offset).coords[0]

    f = 0.5 * width
    line = [(pt[0] + np.cos(angle) * f, pt[1] + np.sin(angle) * f), (pt[0] - np.cos(angle) * f, pt[1] - np.sin(angle) * f)]

    return line

points_in_polygon(points: np.ndarray, polygon: Polygon) -> np.ndarray

Determine points that are inside a polygon, taking holes into account.

Parameters

points : numpy.array Nx2 - array polygon : shapely.geometry.Polygon Polygon (can have holes)

Source code in hydrolib/dhydamo/geometry/spatial.py

def points_in_polygon(points: np.ndarray, polygon: Polygon) -> np.ndarray:
    """
    Determine points that are inside a polygon, taking
    holes into account.

    Parameters
    ----------
    points : numpy.array
        Nx2 - array
    polygon : shapely.geometry.Polygon
        Polygon (can have holes)
    """
    # First select points in square box around polygon
    ptx, pty = points.T
    mainindex = possibly_intersecting(
        dataframebounds=np.c_[[ptx, pty, ptx, pty]], geometry=polygon)
    boxpoints = points[mainindex]

    extp = path.Path(polygon.exterior)
    intps = [path.Path(interior) for interior in polygon.interiors]

    # create first index. Everything within exterior is True
    index = extp.contains_points(boxpoints)

    # set points in holes also to nan
    if intps:
        subset = boxpoints[index]
        # Start with all False
        subindex = np.zeros(len(subset), dtype=bool)

        for intp in intps:
            # update mask, set to True where point in interior
            subindex = subindex | intp.contains_points(subset)

        # Everything within interiors should be True
        # So, set everything within interiors (subindex == True), to True
        index[np.where(index)[0][subindex]] = False

    # Set index in main index to False
    mainindex[np.where(mainindex)[0][~index]] = False

    return mainindex

possibly_intersecting(dataframebounds, geometry, buffer=0)

Finding intersecting profiles for each branch is a slow process in case of large datasets To speed this up, we first determine which profile intersect a square box around the branch With the selection, the interseting profiles can be determines much faster.

Parameters

dataframebounds : numpy.array geometry : shapely.geometry.Polygon

Source code in hydrolib/dhydamo/geometry/spatial.py

def possibly_intersecting(dataframebounds, geometry, buffer=0):
    """
    Finding intersecting profiles for each branch is a slow process in case of large datasets
    To speed this up, we first determine which profile intersect a square box around the branch
    With the selection, the interseting profiles can be determines much faster.

    Parameters
    ----------
    dataframebounds : numpy.array
    geometry : shapely.geometry.Polygon
    """

    geobounds = geometry.bounds
    idx = (
        (dataframebounds[0] - buffer < geobounds[2]) &
        (dataframebounds[2] + buffer > geobounds[0]) &
        (dataframebounds[1] - buffer < geobounds[3]) &
        (dataframebounds[3] + buffer > geobounds[1])
    )
    # Get intersecting profiles
    return idx

rotate_coordinates(origin, theta, xcrds, ycrds)

Rotate coordinates around origin (x0, y0) with a certain angle (radians)

Source code in hydrolib/dhydamo/geometry/spatial.py

def rotate_coordinates(origin, theta, xcrds, ycrds):
    """
    Rotate coordinates around origin (x0, y0) with a certain angle (radians)
    """
    x0, y0 = origin
    xcrds_rot = x0 + (xcrds - x0) * np.cos(theta) + (ycrds - y0) * np.sin(theta)
    ycrds_rot = y0 - (xcrds - x0) * np.sin(theta) + (ycrds - y0) * np.cos(theta)
    return xcrds_rot, ycrds_rot

viz