Clip WBD to 3DEP DEM domain

data/usgs/acquire_and_preprocess_3dep_dems.py

@@ -7,6 +7,7 @@
 import subprocess
 import sys
 import traceback
+import geopandas as gpd
 
 from concurrent.futures import ProcessPoolExecutor, as_completed, wait
 from datetime import datetime
@@ -104,6 +105,8 @@ def acquire_and_preprocess_3dep_dems(extent_file_path,
 
     # download dems, setting projection, block size, etc
     __download_usgs_dems(extent_file_names, target_output_folder_path, number_of_jobs, retry)
+
+    polygonize(target_output_folder_path)
 
     end_time = datetime.now()
     fh.print_end_header('Loading 3dep dems', start_time, end_time)
@@ -273,6 +276,53 @@ def download_usgs_dem_file(extent_file,
         logging.info(msg)
 
 
+def polygonize(target_output_folder_path):
+    """
+    Create a polygon of 3DEP domain from individual HUC6 DEMS which are then dissolved into a single polygon
+    """
+    dem_domain_file = os.path.join(target_output_folder_path, 'HUC6_dem_domain.gpkg')
+
+    msg = f" - Polygonizing -- {dem_domain_file} - Started"
+    print(msg)
+    logging.info(msg)
+
+    dem_files = glob.glob(os.path.join(target_output_folder_path, '*_dem.tif'))
+    dem_gpkgs = gpd.GeoDataFrame()
+
+    for n, dem_file in enumerate(dem_files):
+        edge_tif = f'{os.path.splitext(dem_file)[0]}_edge.tif'
+        edge_gpkg = f'{os.path.splitext(edge_tif)[0]}.gpkg'
+
+        # Calculate a constant valued raster from valid DEM cells
+        if not os.path.exists(edge_tif):
+            subprocess.run(['gdal_calc.py', '-A', dem_file, f'--outfile={edge_tif}', '--calc=where(A > -900, 1, 0)', '--co', 'BIGTIFF=YES', '--co', 'NUM_THREADS=ALL_CPUS', '--co', 'TILED=YES', '--co', 'COMPRESS=LZW', '--co', 'SPARSE_OK=TRUE', '--type=Byte', '--quiet'])
+
+        # Polygonize constant valued raster
+        subprocess.run(['gdal_polygonize.py', '-8', edge_tif, '-q', '-f', 'GPKG', edge_gpkg])
+
+        gdf = gpd.read_file(edge_gpkg)
+
+        if n == 0:
+            dem_gpkgs = gdf
+        else:
+            dem_gpkgs = dem_gpkgs.append(gdf)
+
+        os.remove(edge_tif)
+
+    dem_gpkgs['DN'] = 1
+    dem_dissolved = dem_gpkgs.dissolve(by='DN')
+    dem_dissolved.to_file(dem_domain_file, driver='GPKG')
+
+    if not os.path.exists(dem_domain_file):
+        msg = f" - Polygonizing -- {dem_domain_file} - Failed"
+        print(msg)
+        logging.error(msg)
+    else:
+        msg = f" - Polygonizing -- {dem_domain_file} - Complete"
+        print(msg)
+        logging.info(msg)
+
+
 def __setup_logger(output_folder_path):
 
     start_time = datetime.now()

docs/CHANGELOG.md

@@ -1,6 +1,19 @@
 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.14.1 - 2022-12-03 - [PR #753](https://github.com/NOAA-OWP/inundation-mapping/pull/753)
+
+Creates a polygon of 3DEP DEM domain (to eliminate errors caused by stream networks with no DEM data in areas of HUCs that are outside of the U.S. border) and uses the polygon layer to clip the WBD and stream network (to a buffer inside the WBD).
+
+### Additions
+- `data/usgs/acquire_and_preprocess_3dep_dems.py`: Adds creation of 3DEP domain polygon by polygonizing all HUC6 3DEP DEMs and then dissolving them.
+- `src/gms/run_by_unit.sh`: Adds 3DEP domain polygon .gpkg as input to `src/clip_vectors_to_wbd.py`
+
+### Changes
+- `src/clip_vectors_to_wbd.py`: Clips WBD to 3DEP domain polygon and clips streams to a buffer inside the clipped WBD polygon.
+
+<br/><br/>
+
 ## v4.0.14.0 - 2022-12-20 - [PR #769](https://github.com/NOAA-OWP/inundation-mapping/pull/769)
 
 Masks levee-protected areas from the DEM in branch 0 and in highest two stream order branches.

src/clip_vectors_to_wbd.py

@@ -16,6 +16,7 @@ def subset_vector_layers(subset_nwm_lakes,
                          wbd_buffer_filename,
                          wbd_filename,
                          dem_filename,
+                         dem_domain,
                          nwm_lakes,
                          nwm_catchments,
                          subset_nwm_catchments,
@@ -41,6 +42,10 @@ def subset_vector_layers(subset_nwm_lakes,
     wbd_buffer = wbd.copy()
     wbd_buffer.geometry = wbd.geometry.buffer(wbd_buffer_distance, resolution=32)
 
+    # Erase area outside 3DEP domain
+    dem_domain = gpd.read_file(dem_domain)
+    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)
@@ -162,6 +167,7 @@ def subset_vector_layers(subset_nwm_lakes,
                         required=True)
     parser.add_argument('-g','--wbd-filename', help='HUC boundary', required=True)
     parser.add_argument('-i','--dem-filename', help='DEM filename', required=True)
+    parser.add_argument('-j','--dem-domain', help='DEM domain polygon', required=True)
     parser.add_argument('-l','--nwm-lakes', help='NWM Lakes', required=True)    
     parser.add_argument('-m','--nwm-catchments', help='NWM catchments',
                         required=True)	 

src/gms/run_by_unit.sh

@@ -81,6 +81,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+=" -l $input_nwm_lakes"
 cmd_args+=" -m $input_nwm_catchments"
 cmd_args+=" -n $outputHucDataDir/nwm_catchments_proj_subset.gpkg"