diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md
index 84e549954..5bd79bc5f 100644
--- a/docs/CHANGELOG.md
+++ b/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`
+
+
+
## v4.0.17.4 - 2023-01-06 - [PR#781](https://github.com/NOAA-OWP/inundation-mapping/pull/781)
diff --git a/src/clip_vectors_to_wbd.py b/src/clip_vectors_to_wbd.py
index 60f8639a9..b4cd7d27a 100644
--- a/src/clip_vectors_to_wbd.py
+++ b/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')
diff --git a/src/gms/buffer_stream_branches.py b/src/gms/buffer_stream_branches.py
index cb1426d2c..06313f5ad 100755
--- a/src/gms/buffer_stream_branches.py
+++ b/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)
\ No newline at end of file
diff --git a/src/gms/derive_level_paths.py b/src/gms/derive_level_paths.py
index aaf2eeb19..fd74ddbd3 100755
--- a/src/gms/derive_level_paths.py
+++ b/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):
diff --git a/src/gms/mask_dem.py b/src/gms/mask_dem.py
old mode 100644
new mode 100755
index d8916c2fb..42f0a8feb
--- a/src/gms/mask_dem.py
+++ b/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__':
diff --git a/src/gms/remove_error_branches.py b/src/gms/remove_error_branches.py
old mode 100644
new mode 100755
index 7dfaf26f8..30afbc44b
--- a/src/gms/remove_error_branches.py
+++ b/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
diff --git a/src/gms/run_by_unit.sh b/src/gms/run_by_unit.sh
index 38a70d6ac..67bc2e244 100755
--- a/src/gms/run_by_unit.sh
+++ b/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
diff --git a/src/gms/toDo.md b/src/gms/toDo.md
old mode 100644
new mode 100755
diff --git a/src/split_flows.py b/src/split_flows.py
index 9d83122d0..100be57b1 100755
--- a/src/split_flows.py
+++ b/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)