[1pt] PR: Fix erroneous branch inundation in levee-protected areas

docs/CHANGELOG.md

@@ -1,6 +1,20 @@
 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.4.8.3 - 2024-01-05 - [PR#1059](https://github.com/NOAA-OWP/inundation-mapping/pull/1059)
+
+Fixes erroneous branch inundation in levee-protected areas.
+
+Levees disrupt the natural hydrology and can create large catchments that contain low-lying areas in levee-protected areas that are subject to being inundated in the REM (HAND) grid. However, these low-lying areas are hydrologically disconnected from the stream associated with the catchment and can be erroneously inundated. Branch inundation in levee-protected areas is now confined to the catchment for the levelpath.
+
+### Changes
+
+- `src/`
+    - `delineate_hydros_and_produce_HAND.sh`: Adds input argument for catchments.
+    - `mask_dem.py`: Adds DEM masking for areas of levee-protected areas that are not in the levelpath catchment.
+
+<br/><br/>
+
 ## v4.4.8.2 - 2023-12-12 - [PR#1052](https://github.com/NOAA-OWP/inundation-mapping/pull/1052)
 
 The alpha test for v4.4.8.1 came back with a large degradation in skill and we noticed that the global manning's roughness file was changed in v4.4.7.1 - likely in error.

src/delineate_hydros_and_produce_HAND.sh

@@ -23,6 +23,7 @@ if [ "$mask_leveed_area_toggle" = "True" ] && [ -f $tempHucDataDir/LeveeProtecte
     python3 $srcDir/mask_dem.py \
         -dem $tempCurrentBranchDataDir/dem_meters_$current_branch_id.tif \
         -nld $tempHucDataDir/LeveeProtectedAreas_subset.gpkg \
+        -catchments $z_arg \
         -out $tempCurrentBranchDataDir/dem_meters_$current_branch_id.tif \
         -b $branch_id_attribute \
         -i $current_branch_id \

src/mask_dem.py

@@ -5,6 +5,7 @@
 
 import fiona
 import geopandas as gpd
+import numpy as np
 import pandas as pd
 import rasterio as rio
 from rasterio.mask import mask
@@ -14,6 +15,7 @@ def mask_dem(
     dem_filename: str,
     nld_filename: str,
     levee_id_attribute: str,
+    catchments_filename: str,
     out_dem_filename: str,
     branch_id_attribute: str,
     branch_id: int,
@@ -22,7 +24,8 @@ def mask_dem(
 ):
     """
     Masks levee-protected areas from DEM in branch 0 or if the level path is associated with a levee
-    (determined in src/associate_levelpaths_with_levees.py).
+    (determined in src/associate_levelpaths_with_levees.py). Also masks parts of levee-protected areas
+    through which level paths flow that are not in the level path catchment.
 
     Parameters
     ----------
@@ -41,35 +44,39 @@ def mask_dem(
     branch_zero_id: int
         Branch 0 ID number
     levee_levelpaths: str
-        Path to levee-levelpath association file.
+        Path to levee-levelpath association file (generated by src/associate_levelpaths_with_levees.py)
     """
 
-    # Rasterize if branch zero
-    if branch_id == branch_zero_id:
-        with rio.open(dem_filename) as dem, fiona.open(nld_filename) as leveed:
-            dem_profile = dem.profile.copy()
+    assert os.path.exists(dem_filename), f"DEM file {dem_filename} does not exist"
+    assert os.path.exists(nld_filename), f"NLD file {nld_filename} does not exist"
+    assert os.path.exists(catchments_filename), f"Catchments file {catchments_filename} does not exist"
 
-            geoms = [feature["geometry"] for feature in leveed]
+    dem_masked = None
 
-            # Mask out levee-protected areas from DEM
-            out_dem_masked, _ = mask(dem, geoms, invert=True)
+    with rio.open(dem_filename) as dem:
+        dem_profile = dem.profile.copy()
+        nodata = dem.nodata
 
-            with rio.open(out_dem_filename, "w", **dem_profile, BIGTIFF='YES') as dest:
-                dest.write(out_dem_masked[0, :, :], indexes=1)
+        if branch_id == branch_zero_id:
+            # Mask if branch zero
+            with fiona.open(nld_filename) as leveed:
+                geoms = [feature["geometry"] for feature in leveed]
 
-    elif os.path.exists(levee_levelpaths):
-        levee_levelpaths = pd.read_csv(levee_levelpaths)
+            if len(geoms) > 0:
+                dem_masked, _ = mask(dem, geoms, invert=True)
 
-        # Select levees associated with branch
-        levee_levelpaths = levee_levelpaths[levee_levelpaths[branch_id_attribute] == branch_id]
+        elif os.path.exists(levee_levelpaths):
+            # Mask levee-protected areas protected against level path
+            levee_levelpaths = pd.read_csv(levee_levelpaths)
 
-        # Get levee IDs
-        levelpath_levees = list(levee_levelpaths[levee_id_attribute])
+            # Select levees associated with branch
+            levee_levelpaths = levee_levelpaths[levee_levelpaths[branch_id_attribute] == branch_id]
 
-        if len(levelpath_levees) > 0:
-            with rio.open(dem_filename) as dem:
+            # Get levee IDs
+            levelpath_levees = list(levee_levelpaths[levee_id_attribute])
+
+            if len(levelpath_levees) > 0:
                 leveed = gpd.read_file(nld_filename)
-                dem_profile = dem.profile.copy()
 
                 # Get geometries of levee protected areas associated with levelpath
                 geoms = [
@@ -78,11 +85,40 @@ def mask_dem(
                     if feature[levee_id_attribute] in levelpath_levees
                 ]
 
-                # Mask out levee-protected areas from DEM
-                out_dem_masked, _ = mask(dem, geoms, invert=True)
+                if len(geoms) > 0:
+                    dem_masked, _ = mask(dem, geoms, invert=True)
+
+            # Mask levee-protected areas not protected against level path
+            catchments = gpd.read_file(catchments_filename)
+            leveed = gpd.read_file(nld_filename)
+
+            leveed_area_catchments = gpd.overlay(catchments, leveed, how="union")
+
+            # Select levee catchments not associated with level path
+            levee_catchments_to_mask = leveed_area_catchments.loc[
+                ~leveed_area_catchments[levee_id_attribute].isna() & leveed_area_catchments['ID'].isna(), :
+            ]
+
+            geoms = [feature["geometry"] for i, feature in levee_catchments_to_mask.iterrows()]
+
+            levee_catchments_masked = None
+            if len(geoms) > 0:
+                levee_catchments_masked, _ = mask(dem, geoms, invert=True)
 
+            out_masked = None
+            if dem_masked is None:
+                if levee_catchments_masked is not None:
+                    out_masked = levee_catchments_masked
+
+            else:
+                if levee_catchments_masked is None:
+                    out_masked = dem_masked
+                else:
+                    out_masked = np.where(levee_catchments_masked == nodata, nodata, dem_masked)
+
+            if out_masked is not None:
                 with rio.open(out_dem_filename, "w", **dem_profile, BIGTIFF='YES') as dest:
-                    dest.write(out_dem_masked[0, :, :], indexes=1)
+                    dest.write(out_masked[0, :, :], indexes=1)
 
 
 if __name__ == '__main__':
@@ -91,8 +127,13 @@ def mask_dem(
     parser.add_argument(
         '-nld', '--nld-filename', help='NLD levee-protected areas filename', required=True, type=str
     )
+    parser.add_argument(
+        '-catchments', '--catchments-filename', help='NWM catchments filename', required=True, type=str
+    )
     parser.add_argument('-l', '--levee-id-attribute', help='Levee ID attribute name', required=True, type=str)
-    parser.add_argument('-out', '--out-dem-filename', help='out DEM filename', required=True, type=str)
+    parser.add_argument(
+        '-out', '--out-dem-filename', help='DEM filename to be written', required=True, type=str
+    )
     parser.add_argument(
         '-b', '--branch-id-attribute', help='Branch ID attribute name', required=True, type=str
     )