Source code for watershed_workflow.hydrography

"""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)