v4.0.18.0 Clip WBD and branch buffer polygons to DEM domain (#780)

docs/CHANGELOG.md

@@ -1,6 +1,24 @@
 All notable changes to this project will be documented in this file.
 We follow the [Semantic Versioning 2.0.0](http://semver.org/) format.
 
+## v4.0.18.0 - 2023-01-03 - [PR#780](https://github.com/NOAA-OWP/inundation-mapping/pull/780)
+
+Clips WBD and stream branch buffer polygons to DEM domain.
+
+### Changes
+
+- `src/`
+ - `clip_vectors_to_wbd.py`: Clips WBD polygon to DEM domain
+
+ - `gms/`
+ - `buffer_stream_branches.py`: Clips branch buffer polygons to DEM domain
+ - `derive_level_paths.py`: Stop processing if no branches exist
+ - `mask_dem.py`: Checks if stream file exists before continuing
+ - `remove_error_branches.py`: Checks if error_branches has data before continuing
+ - `run_by_unit.sh`: Adds DEM domain as bash variable and adds it as an argument to calling `clip_vectors_to_wbd.py` and `buffer_stream_branches.py`
+
+<br/><br/>
+
 
 ## v4.0.17.4 - 2023-01-06 - [PR#781](https://github.com/NOAA-OWP/inundation-mapping/pull/781)
 

src/clip_vectors_to_wbd.py

@@ -37,18 +37,20 @@ def subset_vector_layers(subset_nwm_lakes,
  with rio.open(dem_filename) as dem_raster:
  dem_cellsize = max(dem_raster.res)
 
- # Get wbd buffer
- wbd = gpd.read_file(wbd_filename)
- wbd_buffer = wbd.copy()
- wbd_buffer.geometry = wbd.geometry.buffer(wbd_buffer_distance, resolution=32)
-
  # Erase area outside 3DEP domain
+ wbd = gpd.read_file(wbd_filename)
  dem_domain = gpd.read_file(dem_domain)
+ wbd = gpd.clip(wbd, dem_domain)
+ wbd.to_file(wbd_filename, layer='WBDHU8')
+
+ # Get wbd buffer
+ wbd_buffer = wbd.copy()
+ wbd_buffer.geometry = wbd_buffer.geometry.buffer(wbd_buffer_distance, resolution=32)
  wbd_buffer = gpd.clip(wbd_buffer, dem_domain)
 
  # Make the streams buffer smaller than the wbd_buffer so streams don't reach the edge of the DEM
- wbd_streams_buffer = wbd.copy()
- wbd_streams_buffer.geometry = wbd.geometry.buffer(wbd_buffer_distance-3*dem_cellsize, resolution=32)
+ wbd_streams_buffer = wbd_buffer.copy()
+ wbd_streams_buffer.geometry = wbd_streams_buffer.geometry.buffer(-3*dem_cellsize, resolution=32)
 
  great_lakes = gpd.read_file(great_lakes, mask=wbd_buffer).reset_index(drop=True)
 
@@ -57,11 +59,9 @@ def subset_vector_layers(subset_nwm_lakes,
 
  # Clip excess lake area
  great_lakes = gpd.clip(great_lakes, wbd_buffer)
- great_lakes_streams = gpd.clip(great_lakes, wbd_streams_buffer)
 
  # Buffer remaining lake area
  great_lakes.geometry = great_lakes.buffer(lake_buffer_distance)
- great_lakes_streams.geometry = great_lakes_streams.buffer(lake_buffer_distance)
 
  # Removed buffered GL from WBD buffer
  wbd_buffer = gpd.overlay(wbd_buffer, great_lakes, how='difference')

src/gms/buffer_stream_branches.py

@@ -1,6 +1,7 @@
 #!/usr/bin/env python3
 
 import os
+import geopandas as gpd
 from stream_branches import StreamNetwork
 from stream_branches import StreamBranchPolygons
 import argparse
@@ -10,6 +11,7 @@
 
  # parse arguments
  parser = argparse.ArgumentParser(description='Generates branch polygons')
+ parser.add_argument('-a', '--dem-domain', help='DEM domain file', required=False, type=str)
  parser.add_argument('-s','--streams', help='Streams file to branch', required=True)
  parser.add_argument('-i','--branch-id', help='Attribute with branch ids', required=True)
  parser.add_argument('-d','--buffer-distance', help='Distance to buffer branches to create branch polygons', required=True,type=int)
@@ -19,17 +21,21 @@
  # extract to dictionary
  args = vars(parser.parse_args())
 
- streams_file, branch_id_attribute, buffer_distance, stream_polygons_file, verbose = args["streams"], args["branch_id"] , args["buffer_distance"], args["branches"] , args["verbose"]
+ streams_file, branch_id_attribute, buffer_distance, stream_polygons_file, dem_domain, verbose = args["streams"], args["branch_id"] , args["buffer_distance"], args["branches"], args['dem_domain'] , args["verbose"]
 
  if os.path.exists(streams_file):
  # load file
  stream_network = StreamNetwork.from_file( filename=streams_file,branch_id_attribute=branch_id_attribute,
- values_excluded=None,attribute_excluded=None, verbose = verbose)
+ values_excluded=None, attribute_excluded=None, verbose = verbose)
 
  # make stream polygons 
  stream_polys = StreamBranchPolygons.buffer_stream_branches( stream_network,
  buffer_distance=buffer_distance,
  verbose=verbose )
-
- stream_polys.write(stream_polygons_file,verbose=verbose)
 
+ # Clip to DEM domain
+ if os.path.exists(dem_domain):
+ dem_domain = gpd.read_file(dem_domain)
+ stream_polys.geometry = gpd.clip(stream_polys, dem_domain).geometry
+
+ stream_polys.write(stream_polygons_file,verbose=verbose)

src/gms/derive_level_paths.py

@@ -1,5 +1,6 @@
 #!/usr/bin/env python3
 
+import os
 import argparse
 import geopandas as gpd
 import sys
@@ -23,7 +24,11 @@ def Derive_level_paths(in_stream_network, out_stream_network, branch_id_attribut
  if verbose:
  print("Loading stream network ...")
 
- stream_network = StreamNetwork.from_file(filename=in_stream_network)
+ if os.path.exists(in_stream_network):
+ stream_network = StreamNetwork.from_file(filename=in_stream_network)
+ else:
+ print("Sorry, no branches exist and processing can not continue. This could be an empty file.")
+ sys.exit(FIM_exit_codes.GMS_UNIT_NO_BRANCHES.value) # will send a 60 back
 
  # if there are no reaches at this point
  if (len(stream_network) == 0):
@@ -38,8 +43,7 @@ def Derive_level_paths(in_stream_network, out_stream_network, branch_id_attribut
 
  # values_exluded of 1 and 2 mean where are dropping stream orders 1 and 2. We are leaving those
  # for branch zero.
- stream_network = stream_network.exclude_attribute_values(branch_id_attribute="order_",
- values_excluded=[1,2] )
+ stream_network = stream_network.exclude_attribute_values(branch_id_attribute="order_", values_excluded=[1,2] )
 
  # if there are no reaches at this point (due to filtering)
  if (len(stream_network) == 0):

src/gms/mask_dem.py

@@ -1,5 +1,6 @@
 #!/usr/bin/env python3
 
+import os
 import geopandas as gpd
 import fiona
 import rasterio as rio
@@ -12,22 +13,23 @@ def mask_dem(dem_filename, nld_filename, out_dem_filename, stream_layer, branch_
  """
  Masks levee-protected areas from DEM in branch 0 or if stream order is at least the minimum order (max - 1)
  """
- streams_df = gpd.read_file(stream_layer, ignore_geometry=True)
+ if os.path.exists(stream_layer):
+ streams_df = gpd.read_file(stream_layer, ignore_geometry=True)
 
- # Rasterize if branch zero or if stream order is at least the minimum order (max - 1)
- if (branch_id == branch_zero_id) or (streams_df.loc[streams_df[branch_id_attribute].astype(int)==branch_id, order_attribute].max() >= streams_df[order_attribute].max() - 1):
+  # Rasterize if branch zero or if stream order is at least the minimum order (max - 1)
+  if (branch_id == branch_zero_id) or (streams_df.loc[streams_df[branch_id_attribute].astype(int)==branch_id, order_attribute].max() >= streams_df[order_attribute].max() - 1):
 
- with rio.open(dem_filename) as dem, fiona.open(nld_filename) as leveed:
- dem_data = dem.read(1)
- dem_profile = dem.profile.copy()
+  with rio.open(dem_filename) as dem, fiona.open(nld_filename) as leveed:
+  dem_data = dem.read(1)
+  dem_profile = dem.profile.copy()
 
- geoms = [feature["geometry"] for feature in leveed]
+  geoms = [feature["geometry"] for feature in leveed]
 
- # Mask out levee-protected areas from DEM
- out_dem_masked, _ = mask(dem, geoms, invert=True)
+  # Mask out levee-protected areas from DEM
+  out_dem_masked, _ = mask(dem, geoms, invert=True)
 
- with rio.open(out_dem_filename, "w", **dem_profile, BIGTIFF='YES') as dest:
- dest.write(out_dem_masked[0,:,:], indexes=1)
+  with rio.open(out_dem_filename, "w", **dem_profile, BIGTIFF='YES') as dest:
+  dest.write(out_dem_masked[0,:,:], indexes=1)
 
 
 if __name__ == '__main__':

src/gms/remove_error_branches.py

@@ -68,7 +68,7 @@ def remove_error_branches(logfile, gms_inputs):
  error_branches = pd.concat([error_branches, tmp_df])
 
  # Save list of removed branches
- if len(error_branches) > 0:
+ if error_branches is not None and len(error_branches) > 0:
  pd.DataFrame(error_branches).to_csv(gms_inputs_removed, header=False, index=False)
 
  # Overwrite gms_inputs.csv with error branches removed

src/gms/run_by_unit.sh

@@ -48,6 +48,7 @@ input_NHD_WBHD_layer=WBDHU$hucUnitLength
 
 default_projection_crs="ESRI:102039"
 input_DEM=$inputDataDir/3dep_dems/10m_5070/fim_seamless_3dep_dem_10m_5070.vrt
+input_DEM_domain=$inputDataDir/3dep_dems/10m_5070/HUC6_dem_domain.gpkg
 input_NLD=$inputDataDir/nld_vectors/huc2_levee_lines/nld_preprocessed_"$huc2Identifier".gpkg
 input_bathy_bankfull=$inputDataDir/$bankfull_input_table
 
@@ -81,7 +82,7 @@ cmd_args+=" -e $outputHucDataDir/nwm_headwater_points_subset.gpkg"
 cmd_args+=" -f $outputHucDataDir/wbd_buffered.gpkg"
 cmd_args+=" -g $outputHucDataDir/wbd.gpkg"
 cmd_args+=" -i $input_DEM"
-cmd_args+=" -j $inputDataDir/3dep_dems/10m_5070/HUC6_dem_domain.gpkg"
+cmd_args+=" -j $input_DEM_domain"
 cmd_args+=" -l $input_nwm_lakes"
 cmd_args+=" -m $input_nwm_catchments"
 cmd_args+=" -n $outputHucDataDir/nwm_catchments_proj_subset.gpkg"
@@ -102,7 +103,7 @@ Tcount
 python3 $srcDir/clip_vectors_to_wbd.py $cmd_args
 
 : '
-python3 $srcDir/clip_vectors_to_wbd.py -d $hucNumber -w $input_nwm_flows -l $input_nwm_lakes -r $input_NLD -g $outputHucDataDir/wbd.gpkg -f $outputHucDataDir/wbd_buffered.gpkg -m $input_nwm_catchments -y $input_nwm_headwaters -v $input_LANDSEA -lpf $input_nld_levee_protected_areas -z $outputHucDataDir/nld_subset_levees.gpkg -a $outputHucDataDir/nwm_lakes_proj_subset.gpkg -n $outputHucDataDir/nwm_catchments_proj_subset.gpkg -e $outputHucDataDir/nwm_headwater_points_subset.gpkg -b $outputHucDataDir/nwm_subset_streams.gpkg -x $outputHucDataDir/LandSea_subset.gpkg -lps $outputHucDataDir/LeveeProtectedAreas_subset.gpkg -gl $input_GL_boundaries -lb $lakes_buffer_dist_meters -wb $wbd_buffer -i $input_DEM
+python3 $srcDir/clip_vectors_to_wbd.py -d $hucNumber -w $input_nwm_flows -l $input_nwm_lakes -r $input_NLD -g $outputHucDataDir/wbd.gpkg -f $outputHucDataDir/wbd_buffered.gpkg -m $input_nwm_catchments -y $input_nwm_headwaters -v $input_LANDSEA -lpf $input_nld_levee_protected_areas -z $outputHucDataDir/nld_subset_levees.gpkg -a $outputHucDataDir/nwm_lakes_proj_subset.gpkg -n $outputHucDataDir/nwm_catchments_proj_subset.gpkg -e $outputHucDataDir/nwm_headwater_points_subset.gpkg -b $outputHucDataDir/nwm_subset_streams.gpkg -x $outputHucDataDir/LandSea_subset.gpkg -lps $outputHucDataDir/LeveeProtectedAreas_subset.gpkg -gl $input_GL_boundaries -lb $lakes_buffer_dist_meters -wb $wbd_buffer -i $input_DEM -j $input_DEM_domain
 '
 
 ## Clip WBD8 ##
@@ -118,7 +119,6 @@ date -u
 Tstart
 $srcDir/gms/derive_level_paths.py -i $outputHucDataDir/nwm_subset_streams.gpkg -b $branch_id_attribute -r "ID" -o $outputHucDataDir/nwm_subset_streams_levelPaths.gpkg -d $outputHucDataDir/nwm_subset_streams_levelPaths_dissolved.gpkg -e $outputHucDataDir/nwm_headwaters.gpkg -c $outputHucDataDir/nwm_catchments_proj_subset.gpkg -t $outputHucDataDir/nwm_catchments_proj_subset_levelPaths.gpkg -n $outputHucDataDir/nwm_subset_streams_levelPaths_dissolved_headwaters.gpkg -v -w $outputHucDataDir/nwm_lakes_proj_subset.gpkg
 
-
 # test if we received a non-zero code back from derive_level_paths.py
 subscript_exit_code=$?
 # we have to retrow it if it is not a zero (but it will stop further execution in this script)
@@ -129,7 +129,7 @@ Tcount
 echo -e $startDiv"Generating Stream Branch Polygons for $hucNumber"$stopDiv
 date -u
 Tstart
-$srcDir/gms/buffer_stream_branches.py -s $outputHucDataDir/nwm_subset_streams_levelPaths_dissolved.gpkg -i $branch_id_attribute -d $branch_buffer_distance_meters -b $outputHucDataDir/branch_polygons.gpkg -v 
+$srcDir/gms/buffer_stream_branches.py -a $input_DEM_domain -s $outputHucDataDir/nwm_subset_streams_levelPaths_dissolved.gpkg -i $branch_id_attribute -d $branch_buffer_distance_meters -b $outputHucDataDir/branch_polygons.gpkg -v
 Tcount
 
 ## CREATE BRANCHID LIST FILE

src/split_flows.py

@@ -108,9 +108,15 @@ def split_flows(max_length,
  # split at HUC8 boundaries
  print ('splitting stream segments at HUC8 boundaries')
  flows = gpd.overlay(flows, wbd8, how='union').explode().reset_index(drop=True)
+ flows = flows[~flows.is_empty]
+
+ if (len(flows) == 0):
+ # this is not an exception, but a custom exit code that can be trapped
+ print("No relevant streams within HUC boundaries.")
+ sys.exit(FIM_exit_codes.NO_FLOWLINES_EXIST.value) # will send a 61 back
 
  # check for lake features
- if lakes is not None:
+ if lakes is not None and len(flows) > 0 :
  if len(lakes) > 0:
  print ('splitting stream segments at ' + str(len(lakes)) + ' waterbodies')
  #create splits at lake boundaries
@@ -125,6 +131,11 @@ def split_flows(max_length,
  # remove empty geometries
  flows = flows.loc[~flows.is_empty,:]
 
+ if (len(flows) == 0):
+ # this is not an exception, but a custom exit code that can be trapped
+ print("No relevant streams within HUC boundaries.")
+ sys.exit(FIM_exit_codes.NO_FLOWLINES_EXIST.value) # will send a 61 back
+
  for i,lineString in tqdm(enumerate(flows.geometry),total=len(flows.geometry)):
  # Reverse geometry order (necessary for BurnLines)
  lineString = LineString(lineString.coords[::-1])
@@ -266,7 +277,7 @@ def split_flows(max_length,
  max_length = float(environ['max_split_distance_meters'])
  slope_min = float(environ['slope_min'])
  lakes_buffer_input = float(environ['lakes_buffer_dist_meters'])
-
+ 
  # Parse arguments.
  parser = argparse.ArgumentParser(description='splitflows.py')
  parser.add_argument('-f', '--flows-filename', help='flows-filename',required=True)