From 22ba0df00e403681e12d8f247e35da124536a74a Mon Sep 17 00:00:00 2001 From: Brian Avant Date: Fri, 7 May 2021 12:51:19 -0500 Subject: [PATCH] Remove Great Lakes from wbd buffer - The gl_water_polygons.gpkg layer is used to mask out Great Lakes boundaries and remove NHDPlus HR coastline segments. These segments are causing issues later in run_by_unit.sh and unnecessarily increasing total processing time. Resolves issue #374 --- CHANGELOG.md | 8 ++++++++ fim_run.sh | 2 +- src/clip_vectors_to_wbd.py | 37 ++++++++++++++++++++++++++++++++++--- src/run_by_unit.sh | 12 ++---------- 4 files changed, 45 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e10fd0822..e24396a88 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,13 @@ All notable changes to this project will be documented in this file. We follow the [Semantic Versioning 2.0.0](http://semver.org/) format. +

+## v3.0.15.10 - 2021-05-06 - [PR #375](https://github.com/NOAA-OWP/cahaba/pull/375) + +Remove Great Lakes coastlines from WBD buffer. + +## Changes +- `gl_water_polygons.gpkg` layer is used to mask out Great Lakes boundaries and remove NHDPlus HR coastline segments. +

## v3.0.15.9 - 2021-05-03 - [PR #372](https://github.com/NOAA-OWP/cahaba/pull/372) diff --git a/fim_run.sh b/fim_run.sh index 2cfc744e2..cf5de36da 100755 --- a/fim_run.sh +++ b/fim_run.sh @@ -116,7 +116,7 @@ export input_nwm_catchments=$inputDataDir/nwm_hydrofabric/nwm_catchments.gpkg export input_nwm_flows=$inputDataDir/nwm_hydrofabric/nwm_flows.gpkg export input_nhd_flowlines=$inputDataDir/nhdplus_vectors_aggregate/agg_nhd_streams_adj.gpkg export input_nhd_headwaters=$inputDataDir/nhdplus_vectors_aggregate/agg_nhd_headwaters_adj.gpkg - +export input_GL_boundaries=$inputDataDir/landsea/gl_water_polygons.gpkg ## Input handling ## $srcDir/check_huc_inputs.py -u "$hucList" diff --git a/src/clip_vectors_to_wbd.py b/src/clip_vectors_to_wbd.py index f29217e76..9c3d73470 100755 --- a/src/clip_vectors_to_wbd.py +++ b/src/clip_vectors_to_wbd.py @@ -7,15 +7,38 @@ from shapely.geometry import MultiPolygon,Polygon,Point from utils.shared_functions import getDriver -def subset_vector_layers(hucCode,nwm_streams_filename,nhd_streams_filename,nwm_lakes_filename,nld_lines_filename,nwm_catchments_filename,nhd_headwaters_filename,landsea_filename,wbd_filename,wbd_buffer_filename,subset_nhd_streams_filename,subset_nld_lines_filename,subset_nwm_lakes_filename,subset_nwm_catchments_filename,subset_nhd_headwaters_filename,subset_nwm_streams_filename,subset_landsea_filename,extent,dissolveLinks=False): +def subset_vector_layers(hucCode,nwm_streams_filename,nhd_streams_filename,nwm_lakes_filename,nld_lines_filename,nwm_catchments_filename,nhd_headwaters_filename,landsea_filename,wbd_filename,wbd_buffer_filename,subset_nhd_streams_filename,subset_nld_lines_filename,subset_nwm_lakes_filename,subset_nwm_catchments_filename,subset_nhd_headwaters_filename,subset_nwm_streams_filename,subset_landsea_filename,extent,great_lakes_filename,wbd_buffer_distance,lake_buffer_distance,dissolveLinks=False): hucUnitLength = len(str(hucCode)) # Get wbd buffer wbd = gpd.read_file(wbd_filename) - wbd_buffer = gpd.read_file(wbd_buffer_filename) + wbd_buffer = wbd.copy() + wbd_buffer.geometry = wbd.geometry.buffer(wbd_buffer_distance,resolution=32) projection = wbd_buffer.crs + great_lakes = gpd.read_file(great_lakes_filename, mask = wbd_buffer).reset_index(drop=True) + + if not great_lakes.empty: + print("Masking Great Lakes for HUC{} {}".format(hucUnitLength,hucCode),flush=True) + + # Clip excess lake area + great_lakes = gpd.clip(great_lakes, wbd_buffer) + + # Buffer remaining lake area + great_lakes.geometry = great_lakes.buffer(lake_buffer_distance) + + # Removed buffered GL from WBD buffer + wbd_buffer = gpd.overlay(wbd_buffer, great_lakes, how='difference') + wbd_buffer = wbd_buffer[['geometry']] + wbd_buffer.to_file(wbd_buffer_filename,driver=getDriver(wbd_buffer_filename),index=False) + + else: + wbd_buffer = wbd_buffer[['geometry']] + wbd_buffer.to_file(wbd_buffer_filename,driver=getDriver(wbd_buffer_filename),index=False) + + del great_lakes + # Clip ocean water polygon for future masking ocean areas (where applicable) landsea = gpd.read_file(landsea_filename, mask = wbd_buffer) if not landsea.empty: @@ -25,6 +48,7 @@ def subset_vector_layers(hucCode,nwm_streams_filename,nhd_streams_filename,nwm_l # Find intersecting lakes and writeout print("Subsetting NWM Lakes for HUC{} {}".format(hucUnitLength,hucCode),flush=True) nwm_lakes = gpd.read_file(nwm_lakes_filename, mask = wbd_buffer) + nwm_lakes = nwm_lakes.loc[nwm_lakes.Shape_Area < 18990454000.0] if not nwm_lakes.empty: # Perform fill process to remove holes/islands in the NWM lake polygons @@ -59,6 +83,7 @@ def subset_vector_layers(hucCode,nwm_streams_filename,nhd_streams_filename,nwm_l # Subset nhd streams print("Querying NHD Streams for HUC{} {}".format(hucUnitLength,hucCode),flush=True) nhd_streams = gpd.read_file(nhd_streams_filename, mask = wbd_buffer) + if extent == 'MS': nhd_streams = nhd_streams.loc[nhd_streams.mainstem==1] @@ -120,6 +145,9 @@ def subset_vector_layers(hucCode,nwm_streams_filename,nhd_streams_filename,nwm_l parser.add_argument('-b','--subset-nwm-streams',help='NWM streams subset',required=True) parser.add_argument('-x','--subset-landsea',help='LandSea subset',required=True) parser.add_argument('-extent','--extent',help='FIM extent',required=True) + parser.add_argument('-gl','--great-lakes-filename',help='Great Lakes layer',required=True) + parser.add_argument('-wb','--wbd-buffer-distance',help='WBD Mask buffer distance',required=True,type=int) + parser.add_argument('-lb','--lake-buffer-distance',help='Great Lakes Mask buffer distance',required=True,type=int) parser.add_argument('-o','--dissolve-links',help='remove multi-line strings',action="store_true",default=False) @@ -143,6 +171,9 @@ def subset_vector_layers(hucCode,nwm_streams_filename,nhd_streams_filename,nwm_l subset_nwm_streams_filename = args['subset_nwm_streams'] subset_landsea_filename = args['subset_landsea'] extent = args['extent'] + great_lakes_filename = args['great_lakes_filename'] + wbd_buffer_distance = args['wbd_buffer_distance'] + lake_buffer_distance = args['lake_buffer_distance'] dissolveLinks = args['dissolve_links'] - subset_vector_layers(hucCode,nwm_streams_filename,nhd_streams_filename,nwm_lakes_filename,nld_lines_filename,nwm_catchments_filename,nhd_headwaters_filename,landsea_filename,wbd_filename,wbd_buffer_filename,subset_nhd_streams_filename,subset_nld_lines_filename,subset_nwm_lakes_filename,subset_nwm_catchments_filename,subset_nhd_headwaters_filename,subset_nwm_streams_filename,subset_landsea_filename,extent,dissolveLinks) + subset_vector_layers(hucCode,nwm_streams_filename,nhd_streams_filename,nwm_lakes_filename,nld_lines_filename,nwm_catchments_filename,nhd_headwaters_filename,landsea_filename,wbd_filename,wbd_buffer_filename,subset_nhd_streams_filename,subset_nld_lines_filename,subset_nwm_lakes_filename,subset_nwm_catchments_filename,subset_nhd_headwaters_filename,subset_nwm_streams_filename,subset_landsea_filename,extent,great_lakes_filename,wbd_buffer_distance,lake_buffer_distance,dissolveLinks) diff --git a/src/run_by_unit.sh b/src/run_by_unit.sh index 0c5e65cf5..01cc2c98a 100755 --- a/src/run_by_unit.sh +++ b/src/run_by_unit.sh @@ -37,7 +37,7 @@ input_NLD=$inputDataDir/nld_vectors/huc2_levee_lines/nld_preprocessed_"$huc2Iden # Define the landsea water body mask using either Great Lakes or Ocean polygon input # if [[ $huc2Identifier == "04" ]] ; then - input_LANDSEA=$inputDataDir/landsea/gl_water_polygons.gpkg + input_LANDSEA=$input_GL_boundaries echo -e "Using $input_LANDSEA for water body mask (Great Lakes)" else input_LANDSEA=$inputDataDir/landsea/water_polygons_us.gpkg @@ -51,20 +51,12 @@ Tstart ogr2ogr -f GPKG $outputHucDataDir/wbd.gpkg $input_WBD_gdb $input_NHD_WBHD_layer -where "HUC$hucUnitLength='$hucNumber'" Tcount -## BUFFER WBD ## -echo -e $startDiv"Buffer WBD $hucNumber"$stopDiv -date -u -Tstart -[ ! -f $outputHucDataDir/wbd_buffered.gpkg ] && \ -ogr2ogr -f GPKG -dialect sqlite -sql "select ST_buffer(geom, $wbd_buffer) from 'WBDHU$hucUnitLength'" $outputHucDataDir/wbd_buffered.gpkg $outputHucDataDir/wbd.gpkg -Tcount - ## Subset Vector Layers ## echo -e $startDiv"Get Vector Layers and Subset $hucNumber"$stopDiv date -u Tstart [ ! -f $outputHucDataDir/NHDPlusBurnLineEvent_subset.gpkg ] && \ -$srcDir/clip_vectors_to_wbd.py -d $hucNumber -w $input_nwm_flows -s $input_nhd_flowlines -l $input_nwm_lakes -r $input_NLD -g $outputHucDataDir/wbd.gpkg -f $outputHucDataDir/wbd_buffered.gpkg -m $input_nwm_catchments -y $input_nhd_headwaters -v $input_LANDSEA -c $outputHucDataDir/NHDPlusBurnLineEvent_subset.gpkg -z $outputHucDataDir/nld_subset_levees.gpkg -a $outputHucDataDir/nwm_lakes_proj_subset.gpkg -n $outputHucDataDir/nwm_catchments_proj_subset.gpkg -e $outputHucDataDir/nhd_headwater_points_subset.gpkg -b $outputHucDataDir/nwm_subset_streams.gpkg -x $outputHucDataDir/LandSea_subset.gpkg -extent $extent +$srcDir/clip_vectors_to_wbd.py -d $hucNumber -w $input_nwm_flows -s $input_nhd_flowlines -l $input_nwm_lakes -r $input_NLD -g $outputHucDataDir/wbd.gpkg -f $outputHucDataDir/wbd_buffered.gpkg -m $input_nwm_catchments -y $input_nhd_headwaters -v $input_LANDSEA -c $outputHucDataDir/NHDPlusBurnLineEvent_subset.gpkg -z $outputHucDataDir/nld_subset_levees.gpkg -a $outputHucDataDir/nwm_lakes_proj_subset.gpkg -n $outputHucDataDir/nwm_catchments_proj_subset.gpkg -e $outputHucDataDir/nhd_headwater_points_subset.gpkg -b $outputHucDataDir/nwm_subset_streams.gpkg -x $outputHucDataDir/LandSea_subset.gpkg -extent $extent -gl $input_GL_boundaries -lb $lakes_buffer_dist_meters -wb $wbd_buffer Tcount if [ "$extent" = "MS" ]; then