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)