Merge branch 'dev' into dev-acquire-3dep-bug

Dockerfile

@@ -1,5 +1,5 @@
 ## Temporary image to build the libraries and only save the needed artifacts
-FROM ghcr.io/osgeo/gdal:ubuntu-full-3.4.3 AS builder
+FROM ghcr.io/osgeo/gdal:ubuntu-full-3.8.0 AS builder
 WORKDIR /opt/builder
 ARG dataDir=/data
 ARG projectDir=/foss_fim
@@ -15,7 +15,7 @@ RUN git clone https://github.com/dtarb/taudem.git
 RUN git clone https://github.com/fernandoa123/cybergis-toolkit.git taudem_accelerated_flowDirections
 
 RUN apt-get update --fix-missing && apt-get install -y cmake mpich \
-    libgtest-dev libboost-test-dev libnetcdf-dev openjdk-17-jdk && rm -rf /var/lib/apt/lists/*
+    libgtest-dev libboost-test-dev libnetcdf-dev && rm -rf /var/lib/apt/lists/*
 
 ## Compile Main taudem repo ##
 RUN mkdir -p taudem/bin
@@ -36,24 +36,20 @@ RUN mkdir -p $taudemDir
 RUN mkdir -p $taudemDir2
 
 ## Move needed binaries to the next stage of the image
-RUN cd taudem/bin && mv -t $taudemDir flowdircond aread8 threshold streamnet gagewatershed catchhydrogeo dinfdistdown
+RUN cd taudem/bin && mv -t $taudemDir flowdircond streamnet gagewatershed catchhydrogeo dinfdistdown
 RUN cd taudem_accelerated_flowDirections/taudem/build/bin && mv -t $taudemDir2 d8flowdir dinfflowdir
 
 
-
-
 ###############################################################################################
-
-
-
 # Base Image that has GDAL, PROJ, etc
-FROM ghcr.io/osgeo/gdal:ubuntu-full-3.4.3
+FROM ghcr.io/osgeo/gdal:ubuntu-full-3.8.0
 ARG dataDir=/data
 ENV projectDir=/foss_fim
 ARG depDir=/dependencies
 ENV inputsDir=$dataDir/inputs
 ENV outputsDir=/outputs
 ENV srcDir=$projectDir/src
+ENV toolsDir=$projectDir/tools
 ENV workDir=/fim_temp
 ENV taudemDir=$depDir/taudem/bin
 ENV taudemDir2=$depDir/taudem_accelerated_flowDirections/taudem/build/bin
@@ -71,7 +67,7 @@ RUN mkdir -p $depDir
 COPY --from=builder $depDir $depDir
 
 RUN apt update --fix-missing
-RUN DEBIAN_FRONTEND=noninteractive TZ=Etc/UTC apt install -y p7zip-full python3-pip time mpich=3.3.2-2build1 parallel=20161222-1.1 libgeos-dev=3.8.0-1build1 expect=5.45.4-2build1 tmux rsync tzdata openjdk-17-jdk
+RUN DEBIAN_FRONTEND=noninteractive TZ=Etc/UTC apt install -y p7zip-full python3-pip time mpich parallel libgeos-dev expect tmux rsync tzdata
 
 RUN apt auto-remove
 
@@ -106,7 +102,7 @@ RUN pip3 install pipenv==2022.4.8 && PIP_NO_CACHE_DIR=off pipenv install --syste
 # We download and unzip it to the same file folder that pip deployed the whitebox library.
 # Whitebox also attempts to always download a folder called testdata regardless of use.
 # We added an empty folder to fake out whitebox_tools.py so it doesn't try to download the folder
-RUN wbox_path=/usr/local/lib/python3.8/dist-packages/whitebox/ && \
+RUN wbox_path=/usr/local/lib/python3.10/dist-packages/whitebox/ && \
     wget -P $wbox_path https://www.whiteboxgeo.com/WBT_Linux/WhiteboxTools_linux_musl.zip && \
     unzip -o $wbox_path/WhiteboxTools_linux_musl.zip -d $wbox_path && \
     cp $wbox_path/WBT/whitebox_tools $wbox_path && \

Pipfile

@@ -43,6 +43,7 @@ flake8-pyproject = "==1.2.3"
 pre-commit = "==3.3.3"
 isort = "==5.12.0"
 urllib3 = "==1.26.18"
+pyflwdir = "*"
 
 [requires]
-python_version = "3.8"
+python_version = "3.10"

docs/CHANGELOG.md

@@ -24,6 +24,59 @@ Also used the opportunity to fix a couple of other minor issues:
 
 <br/><br/>
 
+## v4.4.9.0 - 2024-01-12 - [PR#1058](https://github.com/NOAA-OWP/inundation-mapping/pull/1058)
+
+Upgrades base Docker image to GDAL v3.8.0. In order to upgrade past GDAL v.3.4.3 (see #1029), TauDEM's `aread8` was replaced with a module from the `pyflwdir` Python package.
+
+### Additions
+
+- `src/accumulate_headwaters.py`: Uses `pyflwdir` to accumulate headwaters and threshold and create stream pixels.
+
+### Changes
+
+- `Dockerfile`: Upgrade GDAL from v.3.4.3 to v.3.8.0; remove JDK 17 and TauDEM `aread8` and `threshold`.
+- `Pipfile` and `Pipfile.lock`: Add `pyflwdir`, `pycryptodomex` and upgrade Python version.
+- `src/delineate_hydros_and_produce_HAND.sh`: Add `src/accumulate_headwaters.py` and remove TauDEM `aread8` and `threshold`
+
+<br/><br/>
+
+## v4.4.8.4 - 2024-01-12 - [PR#1061](https://github.com/NOAA-OWP/inundation-mapping/pull/1061)
+
+Adds a post-processing tool to compare crosswalked (conflated) `feature_id`s between NWM stream network to DEM-derived reaches. The tool is run if the `-x` flag is added to `fim_pipeline.sh`. Results are computed for branch 0 and saved in a summary file in the HUC output folder.
+
+### Additions
+
+- `tools/evaluate_crosswalk.py`: evaluates crosswalk accuracy using two methods:
+    - intersections: the number of intersections between streamlines
+    - network (or tree): compares the feature_ids of the immediate upstream segments
+
+### Changes
+
+- `Dockerfile`: added `toolsDir` environment variable
+- `fim_pipeline.sh`: added `-x` flag to run crosswalk evaluation tool
+- `fim_post_processing.sh`: changed hardcoded `/foss_fim/tools` to `toolsDir` environment variable
+- `fim_pre_processing.sh`: added `evaluateCrosswalk` environment variable
+- `src/`
+    - `add_crosswalk.py`: fix bug
+    - `delineate_hydros_and_produce_HAND.sh`: added a call to `verify_crosswalk.py` if evaluateCrosswalk is True.
+
+<br/><br/>
+
+## 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.
@@ -34,6 +87,7 @@ The alpha test for v4.4.8.1 came back with a large degradation in skill and we n
 
 <br/><br/>
 
+
 ## v4.4.8.1 - 2023-12-08 - [PR#1047](https://github.com/NOAA-OWP/inundation-mapping/pull/1047)
 
 Upgrades JDK to v.17.0.9 in Docker image to address security vulnerabilities.

fim_pipeline.sh

@@ -48,6 +48,7 @@ usage()
                             will be skipped.
       -isaws            : If this param is included, AWS objects will be used where possible
                         - Note: This feature is not yet implemented.
+      -x                : If this param is included, the crosswalk will be evaluated.
 
 
     Running 'fim_pipeline.sh' is a quicker process than running all three scripts independently; however,

fim_post_processing.sh

@@ -210,7 +210,7 @@ echo
 echo -e $startDiv"Combining crosswalk tables"
 # aggregate outputs
 Tstart
-python3 /foss_fim/tools/combine_crosswalk_tables.py \
+python3 $toolsDir/combine_crosswalk_tables.py \
     -d $outputDestDir \
     -o $outputDestDir/crosswalk_table.csv
 Tcount

fim_pre_processing.sh

@@ -101,6 +101,9 @@ in
     -isaws)
         isAWS=1
         ;;
+    -x)
+        evaluateCrosswalk=1
+        ;;
     *) ;;
     esac
     shift
@@ -129,6 +132,7 @@ if [ "$jobBranchLimit" = "" ]; then jobBranchLimit=1; fi
 if [ -z "$overwrite" ]; then overwrite=0; fi
 if [ -z "$skipcal" ]; then skipcal=0; fi
 if [ -z "$isAWS" ]; then isAWS=0; fi
+if [ -z "$evaluateCrosswalk" ]; then evaluateCrosswalk=0; fi
 
 # validate and set defaults for the deny lists
 if [ "$deny_unit_list" = "" ]
@@ -245,6 +249,7 @@ echo "export deny_branch_zero_list=$deny_branch_zero_list" >> $args_file
 echo "export has_deny_branch_zero_override=$has_deny_branch_zero_override" >> $args_file
 echo "export isAWS=$isAWS" >> $args_file
 echo "export skipcal=$skipcal" >> $args_file
+echo "export evaluateCrosswalk=$evaluateCrosswalk" >> $args_file
 
 echo "--- Pre-processing is complete"
 

src/accumulate_headwaters.py

@@ -0,0 +1,102 @@
+#!/usr/bin/env python3
+
+import argparse
+import os
+
+import numpy as np
+import pyflwdir
+import rasterio as rio
+
+
+def accumulate_flow(
+    flow_direction_filename,
+    headwaters_filename,
+    flow_accumulation_filename,
+    stream_pixel_filename,
+    flow_accumulation_threshold,
+):
+    """
+    Accumulate headwaters along the flow direction and threshold accumulations to produce stream pixels.
+
+    Parameters
+    ----------
+    flow_direction_filename : str
+        Flow direction filename
+    headwaters_filename : str
+        Headwaters filename
+    flow_accumulation_filename : str
+        Flow accumulation filename
+    stream_pixel_filename : str
+        Stream pixel filename
+    flow_accumulation_threshold : str
+    """
+
+    assert os.path.isfile(flow_direction_filename), 'Flow direction raster does not exist.'
+
+    # Read the flow direction raster
+    with rio.open(flow_direction_filename) as src:
+        data = src.read(1)
+        nodata = src.nodata
+        profile = src.profile
+        # transform = src.transform
+        # crs = src.crs
+        # latlon = crs.to_epsg() == 4326
+
+    # Convert the TauDEM flow direction raster to a pyflwdir flow direction array
+    temp = data.copy()
+
+    temp[data == 1] = 1
+    temp[data == 2] = 128
+    temp[data == 3] = 64
+    temp[data == 4] = 32
+    temp[data == 5] = 16
+    temp[data == 6] = 8
+    temp[data == 7] = 4
+    temp[data == 8] = 2
+    temp[data == nodata] = 247
+
+    temp = temp.astype(np.uint8)
+
+    flw = pyflwdir.from_array(temp, ftype='d8')
+
+    # Read the flow direction raster
+    with rio.open(headwaters_filename) as src:
+        headwaters = src.read(1)
+        nodata = src.nodata
+
+    flowaccum = flw.accuflux(headwaters, nodata=nodata, direction='up')
+
+    stream = np.where(flowaccum > 0, flow_accumulation_threshold, 0)
+
+    # Write the flow accumulation raster
+    profile.update(dtype=flowaccum.dtype)
+    with rio.open(flow_accumulation_filename, 'w', **profile) as dst, rio.open(
+        stream_pixel_filename, 'w', **profile
+    ) as dst2:
+        dst.write(flowaccum, 1)
+        dst2.write(stream, 1)
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        '-fd', '--flow-direction-filename', help='Flow direction filename', required=True, type=str
+    )
+    parser.add_argument('-wg', '--headwaters-filename', help='Headwaters filename', required=True, type=str)
+    parser.add_argument(
+        '-fa', '--flow-accumulation-filename', help='Flow accumulation filename', required=True, type=str
+    )
+    parser.add_argument(
+        '-stream', '--stream-pixel-filename', help='Stream pixel filename', required=True, type=str
+    )
+    parser.add_argument(
+        '-thresh',
+        '--flow-accumulation-threshold',
+        help='Flow accumulation threshold',
+        required=True,
+        type=float,
+    )
+
+    args = parser.parse_args()
+
+    accumulate_flow(**vars(args))

src/add_crosswalk.py

@@ -258,7 +258,11 @@ def add_crosswalk(
             else:
                 update_id = output_flows.loc[output_flows.HydroID == short_id]['HydroID'].item()
 
-            str_order = output_flows.loc[output_flows.HydroID == short_id]['order_'].item()
+            output_order = output_flows.loc[output_flows.HydroID == short_id]['order_']
+            if len(output_order) == 1:
+                str_order = output_order.item()
+            else:
+                str_order = output_order.max()
             sml_segs = pd.concat(
                 [
                     sml_segs,

src/delineate_hydros_and_produce_HAND.sh

@@ -19,6 +19,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 \
@@ -29,13 +30,18 @@ fi
 
 ## D8 FLOW ACCUMULATIONS ##
 echo -e $startDiv"D8 Flow Accumulations $hucNumber $current_branch_id"
+date -u
+Tstart
 $taudemDir/aread8 -p $tempCurrentBranchDataDir/flowdir_d8_burned_filled_$current_branch_id.tif \
     -ad8 $tempCurrentBranchDataDir/flowaccum_d8_burned_filled_$current_branch_id.tif \
     -wg $tempCurrentBranchDataDir/headwaters_$current_branch_id.tif \
     -nc
+Tcount
 
 # THRESHOLD ACCUMULATIONS ##
 echo -e $startDiv"Threshold Accumulations $hucNumber $current_branch_id"
+date -u
+Tstart
 $taudemDir/threshold -ssa $tempCurrentBranchDataDir/flowaccum_d8_burned_filled_$current_branch_id.tif \
     -src $tempCurrentBranchDataDir/demDerived_streamPixels_$current_branch_id.tif \
     -thresh 1
@@ -229,3 +235,4 @@ python3 $srcDir/add_crosswalk.py \
     -k $tempCurrentBranchDataDir/small_segments_$current_branch_id.csv \
     -e $min_catchment_area \
     -g $min_stream_length
+Tcount