[1pt] PR: dev-bathymetric-adjustment

config/params_template.env

@@ -36,7 +36,7 @@ export levee_id_attribute=SYSTEM_ID
 export healed_hand_hydrocondition=true
 
 #### apply bathymetry adjustment to rating curve ####
-export bathymetry_adjust=True
+export bathymetry_adjust="True"
 
 #### estimating bankfull stage in SRCs ####
 # Toggle to run identify_bankfull routine (True=on; False=off)
@@ -48,6 +48,10 @@ export src_subdiv_toggle="True"
 # text to append to output log and hydrotable file names (use for testing/debugging)
 export vrough_suffix=""
 
+#### applying manningN optimization to SRCs ####
+# Toggle to run optimized roughness src routine (True=on; False=off)
+export src_optz_manningN_toggle="False"
+
 #### SRC calibration variables
 #### apply SRC adjustments using USGS rating curve database ####
 # Toggle to run src adjustment routine with USGS rating data (True=on; False=off)

data/wbd/clip_vectors_to_wbd.py

@@ -2,6 +2,7 @@
 
 import argparse
 import logging
+import os
 import sys
 
 import geopandas as gpd
@@ -21,18 +22,26 @@ def extend_outlet_streams(streams, wbd_buffered, wbd):
  Extend outlet streams to nearest buffered WBD boundary
  """
 
- wbd['geometry'] = wbd.geometry.boundary
- wbd = gpd.GeoDataFrame(data=wbd, geometry='geometry')
-
- wbd_buffered["linegeom"] = wbd_buffered.geometry
-
  # Select only the streams that are outlets
  levelpath_outlets = streams[streams['to'] == 0]
 
- # Select only the streams that don't intersect the WBD boundary line
- levelpath_outlets = levelpath_outlets[~levelpath_outlets.intersects(wbd['geometry'].iloc[0])]
+ levelpath_outlets_columns = [x for x in levelpath_outlets.columns]
+
+ # Select streams that intersect the WBD but not the WBD buffer
+ levelpath_outlets = levelpath_outlets.sjoin(wbd)[levelpath_outlets_columns]
+
+ wbd_boundary = wbd.copy()
+ wbd_boundary['geometry'] = wbd_boundary.geometry.boundary
+ wbd_boundary = gpd.GeoDataFrame(data=wbd_boundary, geometry='geometry')
+
+ wbd_buffered["linegeom"] = wbd_buffered.geometry
+
+ levelpath_outlets = levelpath_outlets[
+ ~levelpath_outlets.intersects(wbd_buffered["linegeom"].boundary.iloc[0])
+ ]
 
  levelpath_outlets['nearest_point'] = None
+ levelpath_outlets['nearest_point_wbd'] = None
  levelpath_outlets['last'] = None
 
  levelpath_outlets = levelpath_outlets.explode(index_parts=False)
@@ -45,22 +54,36 @@ def extend_outlet_streams(streams, wbd_buffered, wbd):
  wbd_buffered['geometry'] = wbd_buffered.geometry.boundary
  wbd_buffered = gpd.GeoDataFrame(data=wbd_buffered, geometry='geometry')
 
+ errors = 0
  for index, row in levelpath_outlets.iterrows():
  levelpath_geom = row['last']
  nearest_point = nearest_points(levelpath_geom, wbd_buffered)
+ nearest_point_wbd = nearest_points(levelpath_geom, wbd_boundary.geometry)
 
  levelpath_outlets.at[index, 'nearest_point'] = nearest_point[1]['geometry'].iloc[0]
+ levelpath_outlets.at[index, 'nearest_point_wbd'] = nearest_point_wbd[1].iloc[0]
 
  levelpath_outlets_nearest_points = levelpath_outlets.at[index, 'nearest_point']
+ levelpath_outlets_nearest_points_wbd = levelpath_outlets.at[index, 'nearest_point_wbd']
+
  if isinstance(levelpath_outlets_nearest_points, pd.Series):
  levelpath_outlets_nearest_points = levelpath_outlets_nearest_points.iloc[-1]
-
- levelpath_outlets.at[index, 'geometry'] = LineString(
- list(row['geometry'].coords) + list([levelpath_outlets_nearest_points.coords[0]])
- )
+ if isinstance(levelpath_outlets_nearest_points_wbd, pd.Series):
+ levelpath_outlets_nearest_points_wbd = levelpath_outlets_nearest_points_wbd.iloc[-1]
+
+ # Extend outlet stream if outlet point is outside of the WBD or nearest snap point is within 100m of the WBD boundary
+ outlet_point = Point(row['geometry'].coords[-1])
+ if (outlet_point.distance(levelpath_outlets_nearest_points_wbd) < 100) or (
+ ~outlet_point.intersects(wbd.geometry)[0]
+ ):
+ levelpath_outlets.at[index, 'geometry'] = LineString(
+ list(row['geometry'].coords) + list([levelpath_outlets_nearest_points.coords[0]])
+ )
+ else:
+ errors += 1
 
  levelpath_outlets = gpd.GeoDataFrame(data=levelpath_outlets, geometry='geometry')
- levelpath_outlets = levelpath_outlets.drop(columns=['last', 'nearest_point'])
+ levelpath_outlets = levelpath_outlets.drop(columns=['last', 'nearest_point', 'nearest_point_wbd'])
 
  # Replace the streams in the original file with the extended streams
  streams = streams[~streams['ID'].isin(levelpath_outlets['ID'])]
@@ -116,6 +139,7 @@ def subset_vector_layers(
  logging.info(f"Clip ocean water polygon for {hucCode}")
  landsea = gpd.read_file(landsea, mask=wbd_buffer, engine="fiona")
  if not landsea.empty:
+ os.makedirs(os.path.dirname(subset_landsea), exist_ok=True)
  # print(f"Create landsea gpkg for {hucCode}", flush=True)
  landsea.to_file(
  subset_landsea, driver=getDriver(subset_landsea), index=False, crs=huc_CRS, engine="fiona"
@@ -356,8 +380,11 @@ def subset_vector_layers(
  '-lps', '--subset-levee-protected-areas', help='Levee-protected areas subset', required=True
  )
  parser.add_argument('-osm', '--osm-bridges', help='Open Street Maps gkpg', required=True)
+ parser.add_argument('-osms', '--subset-osm-bridges', help='Open Street Maps subset', required=True)
+ parser.add_argument('-ak', '--is-alaska', help='If in Alaska', action='store_true')
  parser.add_argument('-crs', '--huc-CRS', help='HUC crs', required=True)
 
  args = vars(parser.parse_args())
 
  subset_vector_layers(**args)
+