"""Functions for manipulating combinations of River and SplitHUCs objects"""
from __future__ import annotations
from typing import Dict, List, Any, Tuple
import math
import copy
import logging
import numpy as np
from matplotlib import pyplot as plt
from scipy.spatial import cKDTree
import itertools
import collections
import shapely.ops
import shapely.geometry
import xarray
import watershed_workflow.config
import watershed_workflow.utils
from watershed_workflow.river_tree import River
from watershed_workflow.split_hucs import SplitHUCs
from watershed_workflow.sources import standard_names as names
#
# findOutlets functions
#
[docs]
def findOutletsByCrossings(hucs: SplitHUCs,
                           river: River,
                           tol: float = 10,
                           debug_plot: bool = False) -> None:
    """For each HUC, find all outlets using a river network's crossing points.
    Parameters
    ----------
    hucs : SplitHUCs
        Split HUCs object to find outlets for.
    river : River
        River network to use for finding crossing points.
    tol : float, optional
        Tolerance in map units for clustering crossings, by default 10.
    debug_plot : bool, optional
        Whether to create debug plots showing outlets, by default False.
    """
    # next determine the outlet, and all boundary edges within x m of that outlet
    # note, we avoid recomputing the polygon geometry because it has a
    # notch removed at the outlet, which will cause an error.
    polygons = list(hucs.df['geometry'])
    poly_crossings = []
    for i_sub, poly in enumerate(polygons):
        my_crossings = []
        for reach in river.preOrder():
            if poly.exterior.intersects(reach.linestring):
                my_crossings.append(poly.exterior.intersection(reach.linestring))
        # cluster my_crossings to make sure that multiple crossings are only counted once
        mccl = []
        for crossing in my_crossings:
            mccl.append([crossing.centroid.xy[0][0], crossing.centroid.xy[1][0]])
        my_crossing_centroids = np.array(mccl)
        clusters, cluster_centroids = watershed_workflow.utils.cluster(my_crossing_centroids, tol)
        poly_crossings.append(cluster_centroids)
    logging.info("Crossings by Polygon:")
    for i, c in enumerate(poly_crossings):
        logging.info(f'  Polygon {i}')
        for p in c:
            logging.info(f'    crossing: {p}')
    # unravel the clusters
    all_crossings_l = [c for p in poly_crossings for c in p]
    all_crossings = np.array(all_crossings_l)
    # cluster crossings that are within tolerance across polygons
    crossings_clusters_indices, crossings_clusters_centroids = \
        
watershed_workflow.utils.cluster(all_crossings, tol)
    # now group cluster ids by polygon and polygon ids by cluster
    poly_cluster_indices = dict()
    cluster_poly_indices = collections.defaultdict(list)
    lcv = 0
    for lcv_poly, pc in enumerate(poly_crossings):
        my_inds = []
        for c in pc:
            my_inds.append(crossings_clusters_indices[lcv])
            lcv += 1
        poly_cluster_indices[lcv_poly] = my_inds
        for ci in my_inds:
            cluster_poly_indices[ci].append(lcv_poly)
    # create a tree, recursively finding all polygons with only
    # one crossing -- this must be an outlet -- then removing it
    # from the list, hopefully leaving a downstream polygon with
    # only one outlet.  This must be done N iterations, where N is
    # the maximal number of polygons crossed from 0th order to
    # maximal order.
    logging.info('Constructing outlet list')
    outlets: Dict[int, int] = dict()
    inlets = collections.defaultdict(list)
    itercount = 0
    done = False
    last_outlet = None
    while not done:
        logging.info(f'Iteration = {itercount}')
        logging.info(f'-----------------')
        new_outlets: Dict[int, int] = dict()
        # look for polygons with only one crossing -- this must be an outlet.
        for pi, clusters2 in poly_cluster_indices.items():
            if len(clusters2) == 1 and pi not in outlets:
                # only one crossing cluster, this is the outlet
                cluster_id = clusters2[0]
                new_outlets[pi] = cluster_id
                cluster_poly_indices[cluster_id].remove(pi)
                logging.info(
                    f' poly outlet {pi} : {cluster_id}, {crossings_clusters_centroids[cluster_id]}')
                last_outlet = cluster_id
                last_outlet_poly = pi
        # look for clusters with only one poly -- this must be an inlet
        to_remove = []
        for ci, polys in cluster_poly_indices.items():
            if len(polys) == 1:
                poly_id = polys[0]
                poly_cluster_indices[poly_id].remove(ci)
                logging.info(f' poly inlet {poly_id} : {ci}, {crossings_clusters_centroids[ci]}')
                to_remove.append(ci)
                inlets[poly_id].append(ci)
        for ci in to_remove:
            cluster_poly_indices.pop(ci)
        if debug_plot and len(new_outlets) > 0:
            fig, ax = plt.subplots(1, 1)
            hucs.plot(color='k', ax=ax)
            river.plot(color='b', ax=ax)
            for pi, ci in outlets.items():
                outlet = crossings_clusters_centroids[ci]
                ax.scatter([outlet[0], ], [outlet[1], ], s=100, c='b', marker='o')
            for pi, ci in new_outlets.items():
                outlet = crossings_clusters_centroids[ci]
                ax.scatter([outlet[0], ], [outlet[1], ], s=100, c='r', marker='o')
            for ci in range(len(crossings_clusters_centroids)):
                if ci not in outlets.values() and ci not in new_outlets.values():
                    crossing = crossings_clusters_centroids[ci]
                    ax.scatter([crossing[0], ], [crossing[1], ], s=100, c='k', marker='o')
            from matplotlib import pyplot as plt
            ax.set_title(f'Outlets after iteration {itercount}')
            plt.show()
        outlets.update(new_outlets)
        itercount += 1
        done = itercount > 50 or len(outlets) == len(polygons) or len(new_outlets) == 0
    if last_outlet is not None:
        logging.info(f'last outlet is {last_outlet} in polygon {last_outlet_poly} '
                     'at {crossings_clusters_centroids[last_outlet]}')
    else:
        logging.info(f'did not find a domain outlet')
    # create the output
    outlet_locs = {}
    inlet_locs = {}
    for pi, ci in outlets.items():
        outlet = crossings_clusters_centroids[ci]
        outlet_locs[pi] = shapely.geometry.Point(outlet[0], outlet[1])
    for pi, cis in inlets.items():
        my_inlet_locs = []
        for ci in cis:
            inlet = crossings_clusters_centroids[ci]
            my_inlet_locs.append(shapely.geometry.Point(inlet[0], inlet[1]))
        inlet_locs[pi] = my_inlet_locs
    last_outlet_p = crossings_clusters_centroids[last_outlet]
    last_outlet_loc = shapely.geometry.Point(last_outlet_p[0], last_outlet_p[1])
    hucs.exterior_outlet = last_outlet_loc
    hucs.df[names.OUTLET] = outlet_locs.values() 
[docs]
def findOutletsByElevation(hucs: SplitHUCs, elev_raster: xarray.Dataset) -> None:
    """Find outlets by the minimum elevation on the boundary.
    Parameters
    ----------
    hucs : SplitHUCs
        Split HUCs object to find outlets for.
    elev_raster : xarray.Dataset
        Elevation raster dataset for determining minimum elevations.
    """
    import watershed_workflow
    def _findOutletsByElevation_helper(polygon, crs, elev_raster):
        mesh_points = np.array(polygon.exterior.coords)
        mesh_points_remapped = watershed_workflow.warp.points(mesh_points, crs, elev_raster.crs)
        elevs = watershed_workflow.utils.valuesFromRaster(mesh_points_remapped, elev_raster)
        i = np.argmin(elevs)
        return shapely.geometry.Point(mesh_points[i])
    hucs.exterior_outlet = _findOutletsByElevation_helper(hucs.exterior, hucs.crs, elev_raster)
    outlets = [
        _findOutletsByElevation_helper(poly, hucs.crs, elev_raster) for poly in hucs.polygons()
    ]
    hucs.df[names.OUTLET] = outlets 
[docs]
def findOutletsByHydroseq(hucs: SplitHUCs, river: River, tol: float = 0.0) -> None:
    """Find outlets using the HydroSequence VAA of NHDPlus.
    Finds the minimum hydroseq reach in each HUC, and intersects that
    with the boundary to find the outlet.
    Parameters
    ----------
    hucs : SplitHUCs
        Split HUCs object to find outlets for.
    river : River
        River network with HydroSequence properties.
    tol : float, optional
        Tolerance for buffering reaches, by default 0.0.
    """
    assert river.isHydroseqConsistent()
    polygons = list(hucs.polygons())
    polygon_outlets = [None for poly in hucs.polygons()]
    # iterate over the reaches, sorted by hydrosequence, looking for
    # the first one that intersects the polygon boundary.
    reaches = sorted(river.preOrder(), key=lambda r: r.properties[names.HYDROSEQ])
    if tol > 0:
        reaches = [r.linestring.buffer(tol) for r in reaches]
    else:
        reaches = [r.linestring for r in reaches]
    first = True
    poly_ids = [(i, poly) for (i, poly) in enumerate(polygons)]
    for lcv, reach in enumerate(reaches):
        try:
            j, (poly_i, poly) = next((j,(i,poly)) for (j,(i,poly)) in enumerate(poly_ids) \
                                     
if poly.intersects(reach))
        except StopIteration:
            continue
        else:
            # find the intersection
            logging.debug(f'hydroseq {lcv} is a match for polygon {poly_i}')
            intersect = poly.exterior.intersection(reach)
            if intersect.is_empty:
                # find the nearest point instead
                intersect = shapely.ops.nearest_points(poly.exterior, reach)[0]
            else:
                intersect = intersect.centroid
            if first:
                hucs.exterior_outlet = intersect
                first = False
            polygon_outlets[poly_i] = intersect
            poly_ids.pop(j)
        if len(poly_ids) == 0:
            break
    hucs.df[names.OUTLET] = polygon_outlets 
[docs]
def snapWaterbodies(waterbodies: List[shapely.geometry.base.BaseGeometry], hucs: SplitHUCs,
                    rivers: List[River], tol: float) -> None:
    """Snap waterbodies to HUCs and river linestrings.
    Attempts to make waterbodies that intersect or nearly intersect
    hucs intersect discretely, in that they share common point(s).
    Parameters
    ----------
    waterbodies : List[shapely.geometry.base.BaseGeometry]
        List of waterbody geometries to snap.
    hucs : SplitHUCs
        Split HUCs object containing boundary linestrings.
    rivers : List[River]
        List of river networks to snap to.
    tol : float
        Snapping tolerance in map units.
    """
    # note, this is a fairly fragile algorithm that looks only at
    # discrete points.  A more robust one could be developed if
    # needed.
    mls = shapely.geometry.MultiLineString([r.linestring for river in rivers for r in river]
                                           + [polygon.exterior for polygon in hucs.polygons()])
    for i, wb in enumerate(waterbodies):
        waterbodies[i] = shapely.snap(wb, mls, tol) 
[docs]
def cutAndSnapCrossings(hucs: SplitHUCs, rivers: List[River], tol: float) -> None:
    """Aligns river and HUC objects.
    1. where a reach crosses an external boundary, cut in two and keep
       only internal portion.
    2. where a reach crosses an internal boundary, either:
       - snap the internal boundary to the reach endpoint
       - cut the reach in two
    At the end of the day, we ensure that any place where a river
    crosses a HUC boundary is a discrete point at a both a reach
    endpoint and a HUC boundary segment.  This ensures that those
    points will never move and will always be coincident.
    """
    # check exterior crossings on trunk and leaf
    for river in rivers:
        _cutAndSnapExteriorCrossing(hucs, river, tol)
        for leaf in river.leaf_nodes:
            _cutAndSnapExteriorCrossing(hucs, leaf, tol)
    # check interior crossings on all
    for river in rivers:
        # make a copy as river is modified in-place
        reaches = list(river)
        for reach in reaches:
            _cutAndSnapInteriorCrossing(hucs, reach, tol) 
def _cutAndSnapExteriorCrossing(hucs: SplitHUCs, reach: River, merge_tol: float) -> None:
    """Cut and snap reach at exterior HUC boundary crossings.
    If any reach crosses a HUC boundary:
    1. If it crosses an external boundary, cut the reach in two and
       discard the portion outside of the domain.
    2. If it crosses an internal boundary, ensure there is a
       coordinate of the reach that is on the internal boundary.
    Either way, also ensure there is a coordinate on the HUC
    boundary at the crossing point.
    Modifies no geometry, only topology.
    Parameters
    ----------
    hucs : SplitHUCs
        Split HUCs object containing boundary information.
    reach : River
        River reach to process for boundary crossings.
    merge_tol : float
        Tolerance for merging operations.
    """
    r = reach.linestring
    # first deal with crossings of the HUC exterior boundary -- in
    # this case, the reach linestring gets split in two and the external
    # one is removed.
    for b, spine in hucs.boundaries.items():
        # make a copy as spine is modified in-place
        spine_list = list(spine.items())
        for s, ls_handle in spine_list:
            ls = hucs.linestrings[ls_handle]
            if watershed_workflow.utils.intersects(ls, r):
                logging.info('intersection found')
                new_spine, new_reach = watershed_workflow.utils.cut(ls, r)
                assert len(new_reach) == 1 or len(new_reach) == 2
                assert len(new_spine) == 1 or len(new_spine) == 2
                logging.info("  - cutting reach at external boundary of HUCs:")
                logging.info(f"      split HUC boundary ls into {len(new_spine)} pieces")
                logging.info(f"      split reach ls into {len(new_reach)} pieces")
                # which piece of the reach are we keeping?
                if hucs.exterior.buffer(-1).contains(shapely.geometry.Point(
                        new_reach[0].coords[0])):
                    # keep the upstream (or only) reach ls
                    if len(new_reach) == 2:
                        # confirm other/downstream reach is outside
                        assert not hucs.exterior.contains(
                            shapely.geometry.Point(new_reach[1].coords[-1]))
                    reach.linestring = new_reach[0]
                elif len(new_reach) == 2:
                    if hucs.exterior.buffer(-1).contains(
                            shapely.geometry.Point(new_reach[1].coords[-1])):
                        # keep the downstream reach ls, confirm upstream is outside
                        assert not hucs.exterior.contains(
                            shapely.geometry.Point(new_reach[0].coords[0]))
                        reach.linestring = new_reach[1]
                # keep both pieces of a split huc boundary linestring
                # -- rename the first
                hucs.linestrings[ls_handle] = new_spine[0]
                if len(new_spine) > 1:
                    # -- add the second
                    new_handle = hucs.linestrings.append(new_spine[1])
                    spine.append(new_handle)
def _cutAndSnapInteriorCrossing(hucs: SplitHUCs, reach: River, merge_tol: float) -> None:
    """Cut and snap reach at interior HUC boundary crossings.
    Helper function for cutAndSnapCrossings().
    Modifies no geometry, only topology.
    Parameters
    ----------
    hucs : SplitHUCs
        Split HUCs object containing interior boundary information.
    reach : River
        River reach to process for interior crossings.
    merge_tol : float
        Tolerance for merging operations.
    """
    r = reach.linestring
    # now deal with crossings of the HUC interior boundary -- in this
    # case, the reach linestring cut, then potentially merged to neighbors
    for i, spine in hucs.intersections.items():
        # make a copy as spine is modified in-place
        spine_list = list(spine.items())
        for s, ls_handle in spine_list:
            ls = hucs.linestrings[ls_handle]
            if watershed_workflow.utils.intersects(ls, r):
                new_spine, new_reach = watershed_workflow.utils.cut(ls, r)
                assert len(new_reach) == 1 or len(new_reach) == 2
                assert len(new_spine) == 1 or len(new_spine) == 2
                logging.info("  - snapping reach at internal boundary of HUCs")
                if len(new_reach) == 1:
                    reach.linestring = new_reach[0]
                elif len(new_reach) == 2:
                    logging.info('  branch2')
                    reach.linestring = shapely.geometry.LineString(
                        list(new_reach[0].coords) + list(new_reach[1].coords)[1:])
                    logging.info(f'  seg1: {list(new_reach[0].coords)}')
                    logging.info(f'  seg2: {list(new_reach[1].coords[1:])}')
                    logging.info(
                        f'  splitting at coord: {len(new_reach[0].coords)-1} of {len(reach.linestring.coords)}'
                    )
                    us, ds = reach.split(len(new_reach[0].coords) - 1)
                hucs.linestrings[ls_handle] = new_spine[0]
                if len(new_spine) > 1:
                    assert (len(new_spine) == 2)
                    new_handle = hucs.linestrings.append(new_spine[1])
                    spine.append(new_handle)
[docs]
def snapHUCsJunctions(hucs: SplitHUCs, rivers: List[River], tol: float) -> None:
    """Snaps the junctions of HUC linestrings to endpoints of rivers.
    Modifies HUCs geometry.
    Parameters
    ----------
    hucs : SplitHUCs
        Split HUCs object to modify.
    rivers : List[River]
        List of river networks to snap to.
    tol : float
        Snapping tolerance in map units.
    """
    # make the kdTree of endpoints of all reaches
    coords1 = np.array(
        [reach.linestring.coords[-1] for river in rivers for reach in river.preOrder()])
    coords2 = np.array(
        [reach.linestring.coords[0] for river in rivers for reach in river.leaf_nodes])
    coords = np.concatenate([coords1, coords2], axis=0)
    # limit to x,y
    if (coords.shape[1] != 2):
        coords = coords[:, 0:2]
    kdtree = cKDTree(coords)
    # for each linestring of the HUC spine, find the river outlet that is
    # closest.  If within tolerance, move it
    for ls_handle, ls in hucs.linestrings.items():
        # check point 0, -1
        endpoints = np.array([ls.coords[0], ls.coords[-1]])
        # delete z coord
        if (endpoints.shape[1] != 2):
            endpoints = endpoints[:, 0:2]
        dists, inds = kdtree.query(endpoints)
        if dists.min() < tol:
            new_ls = list(ls.coords)
            if dists[0] < tol:
                new_ls[0] = coords[inds[0]]
                logging.debug(
                    f"  Moving HUC linestring point 0,1: {list(ls.coords)[0]}, {list(ls.coords)[-1]}"
                )
                logging.debug("        point 0 to river at %r" % list(new_ls[0]))
                # zap points until we are outside of the ball -- otherwise we can get weird crossings
                count = 0
                while watershed_workflow.utils.computeDistance(new_ls[1], new_ls[0]) < tol:
                    new_ls.pop(1)
                    count += 1
                logging.debug(f"    also removed {count} other coords")
            if dists[1] < tol:
                new_ls[-1] = coords[inds[1]]
                logging.debug(
                    f"  Moving HUC linestring point 0,1: {list(ls.coords)[0]}, {list(ls.coords)[-1]}"
                )
                logging.debug("        point -1 to river at %r" % list(new_ls[-1]))
                # zap points until we are outside of the ball -- otherwise we can get weird crossings
                count = 0
                while watershed_workflow.utils.computeDistance(new_ls[-2], new_ls[-1]) < tol:
                    new_ls.pop(-2)
                    count += 1
                logging.debug(f"    also removed {count} other coords")
            hucs.linestrings[ls_handle] = shapely.geometry.LineString(new_ls) 
def _findContainingPolygon(hucs: SplitHUCs, linestring: shapely.geometry.LineString) -> int:
    """Find the polygon that contains the largest portion of linestring.
    Parameters
    ----------
    hucs : SplitHUCs
        Split HUCs object containing polygons.
    linestring : shapely.geometry.LineString
        Linestring to find containing polygon for.
    Returns
    -------
    int
        Index of the polygon containing the largest portion of the linestring.
    """
    return max((i for i in range(len(hucs))),
               key=lambda i: linestring.intersection(hucs.computePolygon(i)).length)
[docs]
def snapReachEndpoints(hucs: SplitHUCs, river: River, tol: float) -> None:
    """Snap river endpoints to HUC linestrings and insert that point into the boundary.
    Note this is O(n^2), and could be made more efficient.
    Modifies reach geometry.
    Parameters
    ----------
    hucs : SplitHUCs
        Split HUCs object containing boundary linestrings.
    river : River
        River network to snap endpoints for.
    tol : float
        Snapping tolerance in map units.
    """
    to_add = []
    reaches = list(river)
    # check the river oulet
    reach_ls = river.linestring
    done = False
    for b, component in itertools.chain(hucs.boundaries.items(), hucs.intersections.items()):
        for s, huc_ls_handle in component.items():
            huc_ls = hucs.linestrings[huc_ls_handle]
            logging.debug("  - checking reach coord: %r" % list(reach_ls.coords[-1]))
            logging.debug("  - huc_ls coords: {0}".format(list(huc_ls.coords)))
            # find the nearest point to the endpoint if it is within tol
            new_coord = watershed_workflow.utils.findNearestPoint(reach_ls.coords[-1], huc_ls, tol)
            if new_coord is not None:
                if any(watershed_workflow.utils.isClose(new_coord, c) for c in huc_ls.coords):
                    # the endpoint is already discretely in the HUC
                    # boundary, likely done in a junction snap
                    done = True
                    break
                logging.debug("  - new coord: {0}".format(new_coord))
                logging.debug("  - snapped reach: %r to %r" % (reach_ls.coords[-1], new_coord))
                # keep a list of all points to add, which are all added at once
                to_add.append((huc_ls_handle, component, -1, river))
                # remove points on the reach that are
                # closer to the huc -- this deals with the
                # case that multiple discrete points are
                # on the "wrong" side of the internal
                # boundary.
                coords = list(reach_ls.coords)
                while len(coords) > 2 and \
                      
watershed_workflow.utils.computeDistance(new_coord, coords[-2]) \
                      
< watershed_workflow.utils.computeDistance(new_coord, coords[-1]):
                    coords.pop(-1)
                coords[-1] = new_coord
                reach_ls = shapely.geometry.LineString(coords)
                river.linestring = reach_ls
                # if we add a point on this huc linestring, don't
                # add it to any other
                done = True
                break
        if done: break  # break out of both for loops
    # now check the upstream end of all reaches
    for reach in reaches:
        reach_ls = reach.linestring
        done = False
        # only consider endpoints for whom the touching reaches span multiple polygons
        reach_linestrings = [reach_ls, ] + [
            watershed_workflow.utils.reverseLineString(c.linestring) for c in reach.children
        ]
        touching_polygons = set(_findContainingPolygon(hucs, r) for r in reach_linestrings)
        if len(touching_polygons) > 1:
            # find the component it touches
            for b, component in itertools.chain(hucs.boundaries.items(),
                                                hucs.intersections.items()):
                for s, huc_ls_handle in component.items():
                    huc_ls = hucs.linestrings[huc_ls_handle]
                    logging.debug("  - checking reach coord: %r" % list(reach_ls.coords[0]))
                    logging.debug("  - huc_ls coords: {0}".format(list(huc_ls.coords)))
                    # only consider endpoints for whom one reach
                    # intersects the HUC boundary.  If the endpoint is
                    # close but no reaches intersect the boundary, it
                    # is likely that the endpoint is fully contained
                    # in the polygon.  We shrink the linestring at the
                    # opposite endpoint to avoid counting an
                    # intersection at the other end of the reach.
                    def _shrinkLS(ls):
                        coords = [c for c in ls.coords]
                        coords[-1] = watershed_workflow.utils.computeMidpoint(
                            coords[-2], coords[-1])
                        return shapely.geometry.LineString(coords)
                    #
                    # this is the upstream point of this reach, so
                    # consider this and all children
                    if any(
                            watershed_workflow.utils.intersects(huc_ls, _shrinkLS(ls))
                            for ls in reach_linestrings):
                        # find the nearest point to the endpoint if it is within tol
                        new_coord = watershed_workflow.utils.findNearestPoint(
                            reach_ls.coords[0], huc_ls, tol)
                        if new_coord is not None:
                            if any(
                                    watershed_workflow.utils.isClose(new_coord, c)
                                    for c in huc_ls.coords):
                                # the endpoint is already discretely in the
                                # HUC boundary, likely done as a junction
                                done = True
                                break
                            logging.debug("  - new coord: {0}".format(new_coord))
                            logging.debug(f"    snapped reach: {reach_ls.coords[0]} to {new_coord}")
                            # keep a list of all points to add, which are all added at once
                            to_add.append((huc_ls_handle, component, 0, reach))
                            for ls, r in zip(reach_linestrings, [reach, ] + list(reach.children)):
                                # remove points on the reach that are
                                # closer to the huc -- this deals with the
                                # case that multiple discrete points are
                                # on the "wrong" side of the internal
                                # boundary.
                                coords = list(ls.coords)
                                while len(coords) > 2 and \
                                      
watershed_workflow.utils.computeDistance(new_coord, coords[1]) \
                                      
< watershed_workflow.utils.computeDistance(new_coord, coords[0]):
                                    coords.pop(0)
                                coords[0] = new_coord
                                ls = shapely.geometry.LineString(coords)
                                if r is reach:
                                    r.linestring = ls
                                else:
                                    r.linestring = watershed_workflow.utils.reverseLineString(ls)
                            assert reach.isLocallyContinuous()
                            # if we add a point on this huc linestring, don't
                            # add it to any other
                            done = True
                            break
                if done: break  # break out of two for loops
        if reach.linestring.length == 0.:
            # snapped both endpoints to the same point on the internal
            # boundary, remove the reach note we can safely merge this
            # with either parent or child
            if len(list(reach.siblings)) == 0:
                reach.merge()
            else:
                assert len(reach.children) == 1
                reach.children[0].merge()
    # find the list of points to add to a given linestring
    to_add_dict: Dict[int, List[Any]] = dict()
    for huc_ls_handle, component, endpoint, reach in to_add:
        if huc_ls_handle not in to_add_dict.keys():
            to_add_dict[huc_ls_handle] = list()
        to_add_dict[huc_ls_handle].append((component, endpoint, reach))
    # find the set of points to add to each given linestring
    def isEqual(p1, p2):
        if watershed_workflow.utils.isClose(p1[2].linestring.coords[p1[1]],
                                            p2[2].linestring.coords[p2[1]], 1.e-5):
            assert (p1[0] == p2[0])
            return True
        else:
            return False
    to_add_dict2 = dict()
    for huc_ls_handle, insert_list in to_add_dict.items():
        new_list: List[Any] = []
        for p1 in insert_list:
            if (all(not isEqual(p1, p2) for p2 in new_list)):
                new_list.append(p1)
        to_add_dict2[huc_ls_handle] = new_list
    # add these points to the linestring
    for huc_ls_handle, insert_list in to_add_dict2.items():
        ls = hucs.linestrings[huc_ls_handle]
        # make a list of the coords and a flag to indicate a new
        # coord, then sort it by arclength along the linestring.
        #
        # Note this needs special care if the ls is a loop, or else the endpoint gets sorted twice
        if not watershed_workflow.utils.isClose(ls.coords[0], ls.coords[-1]):
            new_coords = [[p[2].linestring.coords[p[1]], 1] for p in insert_list]
            old_coords = [
                [c, 0] for c in ls.coords
                if not any(watershed_workflow.utils.isClose(c, nc[0], tol) for nc in new_coords)
            ]
            new_ls_coords = sorted(new_coords + old_coords,
                                   key=lambda a: ls.project(shapely.geometry.Point(a[0])))
            # determine the new coordinate indices
            breakpoint_inds = [i for i, (c, f) in enumerate(new_ls_coords) if f == 1]
        else:
            new_coords = [[p[2].linestring.coords[p[1]], 1] for p in insert_list]
            old_coords = [
                [c, 0] for c in ls.coords[:-1]
                if not any(watershed_workflow.utils.isClose(c, nc[0], tol) for nc in new_coords)
            ]
            new_ls_coords = sorted(new_coords + old_coords,
                                   key=lambda a: ls.project(shapely.geometry.Point(a[0])))
            breakpoint_inds = [i for i, (c, f) in enumerate(new_ls_coords) if f == 1]
            assert (len(breakpoint_inds) > 0)
            new_ls_coords = new_ls_coords[breakpoint_inds[0]:] + new_ls_coords[0:breakpoint_inds[0]
                                                                               + 1]
            new_ls_coords[0][1] = 0
            new_ls_coords[-1][1] = 0
            breakpoint_inds = [i for i, (c, f) in enumerate(new_ls_coords) if f == 1]
        # now break into new linestrings
        new_lss = []
        ind_start = 0
        for ind_end in breakpoint_inds:
            assert (ind_end != 0)
            new_lss.append(
                shapely.geometry.LineString([c for (c, f) in new_ls_coords[ind_start:ind_end + 1]]))
            ind_start = ind_end
        assert (ind_start < len(new_ls_coords) - 1)
        new_lss.append(
            shapely.geometry.LineString([tuple(c) for (c, f) in new_ls_coords[ind_start:]]))
        # put all new_lss into the huc list.  Note insert_list[0][0] is the component
        watershed_workflow.utils.logMinMaxMedianSegment(new_lss, "new_lss")
        hucs.linestrings[huc_ls_handle] = new_lss.pop(0)
        new_handles = hucs.linestrings.extend(new_lss)
        insert_list[0][0].extend(new_handles)