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
NOAA-OWP · May 7, 2021 · 22ba0df · 22ba0df
1 parent 25c3ec1
commit 22ba0df
Show file tree

Hide file tree

Showing 4 changed files with 45 additions and 14 deletions.
diff --git 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.
+<br/><br/>
+## 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.
+
 <br/><br/>
 ## 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
@@ -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
@@ -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
@@ -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