From 0466ec4a3cb9f2a24671a2ba3e9ac4087cc2c19f Mon Sep 17 00:00:00 2001 From: MattLuck-NOAA Date: Fri, 6 Jan 2023 14:44:51 -0700 Subject: [PATCH] v4.0.18.0 Clip WBD and branch buffer polygons to DEM domain (#780) --- docs/CHANGELOG.md | 18 ++++++++++++++++++ src/clip_vectors_to_wbd.py | 18 +++++++++--------- src/gms/buffer_stream_branches.py | 14 ++++++++++---- src/gms/derive_level_paths.py | 10 +++++++--- src/gms/mask_dem.py | 24 +++++++++++++----------- src/gms/remove_error_branches.py | 2 +- src/gms/run_by_unit.sh | 8 ++++---- src/gms/toDo.md | 0 src/split_flows.py | 15 +++++++++++++-- 9 files changed, 75 insertions(+), 34 deletions(-) mode change 100644 => 100755 src/gms/mask_dem.py mode change 100644 => 100755 src/gms/remove_error_branches.py mode change 100644 => 100755 src/gms/toDo.md 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)