[8pt] PR: Fix Alaska HydroIDs

docs/CHANGELOG.md

@@ -1,6 +1,24 @@
 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.8.14.0 - 2025-10-10 - [PR#1650](https://github.com/NOAA-OWP/inundation-mapping/pull/1650)
+
+Fixes Alaska HUCs attempt to convert to Int16 if there are too many Hydro IDs and creates error handling for when this issue may arise.
+
+### Changes
+- src/delineate_hydros_and_produce_HAND.sh: Add check for Alaska HUC2 identifier before running conversion to int16 data type.
+- src/process_branch.sh: Add handling for int16 conversion failure error code.
+- src/src_adjust_spatial_obs.py: Add check for Alaska HUC2 identifier to see if process should run in Int16 or Int32.
+- src/utils/fim_enums.py: Add error code for int16 conversion error.
+- tools/catfim/generate_categorical_fim_mapping.py: Process in either in Int16 for rem and HydroIDs, or int32 and float32 respectively in the case of an Alaska HUC.
+- tools/inundation.py: Process in either in Int16 for rem and HydroIDs, or int32 and float32 respectively in the case of an Alaska HUC.
+- tools/mosaic_inundation.py: Add rasterio merge and remove calls to removed script overlapping_inundation.py.
+- tools/convert_to_int16.py: Add raise error for scenarios where catchment raster has too many HydroIDs or HydroIDs that have more than 8 digits.
+
+### Removals
+- /tools/overlapping_inundation.py
+<br />
+
 ## v4.8.13.0 - 2025-10-10 - [PR#1673](https://github.com/NOAA-OWP/inundation-mapping/pull/1673)
 
 This script processes river segment slope data to identify and correct unrealistic values based on user-defined thresholds. It is designed to check against the NWM hydrofab to and fill/extrapolate missing or invalid slope values using valid upstream/downstream values.

src/delineate_hydros_and_produce_HAND.sh

@@ -2,6 +2,7 @@
 
 ## Level is equal to the parent script: 'unit' or 'branch'
 level=$1
+huc2Identifier=${hucNumber:0:2}
 
 if [ "$level" = "branch" ]; then
     b_arg=$tempCurrentBranchDataDir/nwm_subset_streams_levelPaths_$current_branch_id.gpkg
@@ -311,8 +312,12 @@ if [ "$current_branch_id" = "$branch_zero_id" ] && [ "$evaluateCrosswalk" = "1"
         -z $current_branch_id
 fi
 
-## CONVERSION TO INT16 ##
-echo -e $startDiv"Convert GW Catchments and REM to Int16 $hucNumber $current_branch_id"
-python3 $toolsDir/convert_to_int16.py \
-    -b $tempCurrentBranchDataDir
+if [ $huc2Identifier -eq 19 ]; then
+    echo -e "Skipping Int16 Conversion for Alaska HUC"
+else
+    ## CONVERSION TO INT16 ##
+    echo -e $startDiv"Convert GW Catchments and REM to Int16 $hucNumber $current_branch_id"
+    python3 $toolsDir/convert_to_int16.py \
+        -b $tempCurrentBranchDataDir
+fi
 

src/process_branch.sh

@@ -43,6 +43,11 @@ do
         err_exists=1
         echo "***** Branch has no crosswalks *****"
         rm -rf $tempHucDataDir/branches/$branchId/
+    elif [ $code -eq 65]; then
+        echo
+        err_exists=1
+        echo "***** Too many HydroIDs or a HydroID with more than 8 digits in gw catchments to convert to Int16 *****"
+        rm -rf $tempHucDataDir/branches/$branchId/
     elif [ $code -ne 0 ]; then
         echo
         err_exists=1

src/src_adjust_spatial_obs.py

@@ -86,18 +86,31 @@ def process_points(args):
 
     water_edge_df = water_edge_df.to_crs(DEFAULT_FIM_PROJECTION_CRS)
 
-    with open(hydroid_prefixpath, 'r') as file:
-        hydroid_prefix = file.read()
-        int_hid_prefix = int(hydroid_prefix) * 10000
+    process_int16 = True
+    try:
+        with open(hydroid_prefixpath, 'r') as file:
+            hydroid_prefix = file.read()
+            int_hid_prefix = int(hydroid_prefix) * 10000
+    except FileNotFoundError:
+        process_int16 = False
+        int_hid_prefix = 0
 
     ## Use point geometry to determine HAND raster pixel values.
     with rasterio.open(hand_path) as hand_src, rasterio.open(catchments_path) as catchments_src:
-        water_edge_df['hand'] = [np.float32(h[0]) / 1000 for h in hand_src.sample(coords)]
+
+        water_edge_df['hand'] = (
+            [np.float32(h[0]) / 1000 for h in hand_src.sample(coords)]
+            if process_int16
+            else [h[0] for h in hand_src.sample(coords)]
+        )
+
         hydroids = []
 
         for c in catchments_src.sample(coords):
-            hid = int_hid_prefix * -1 + c[0] if c[0] < 0 else int_hid_prefix + c[0]
+            c[0] = c[0] if c[0] >= 0 else c[0] * -1
+            hid = int_hid_prefix + c[0] if process_int16 else c[0]
             hydroids.append(hid)
+
         water_edge_df['hydroid'] = hydroids
 
     water_edge_df = water_edge_df[

src/utils/fim_enums.py

@@ -25,8 +25,13 @@ class FIM_exit_codes(Enum):
                              https://docs.python.org/3/library/enum.html
     '''
 
+    # Codes beginning with 6 are acceptable errors
     UNIT_NO_BRANCHES = 60
     NO_FLOWLINES_EXIST = 61
     EXCESS_UNIT_ERRORS = 62
     NO_BRANCH_LEVELPATHS_EXIST = 63
     NO_VALID_CROSSWALKS = 64
+    CANNOT_CONVERT_HYDROIDS_TO_INT16 = 65
+
+    # @TODO Make system codes that shut down the entire HUC process
+    # Codes beginning with 5 are not acceptable errors and shutdown HUC processing

tools/catfim/generate_categorical_fim_mapping.py

@@ -73,7 +73,9 @@ def produce_stage_based_lid_tifs(
 
     # Subtract HAND gage elevation from HAND WSE to get HAND stage.
     hand_stage_m = datum_adj_wse_m - lid_usgs_elev
-    hand_stage = round(hand_stage_m * 1000)  # convert to mm to match HAND
+    hand_stage = (
+        hand_stage_m if str(huc)[:2] == '19' else round(hand_stage_m * 1000)
+    )  # convert to mm to match HAND
 
     # If hand_stage is negative, write message and exit out
     if hand_stage < 0:
@@ -324,15 +326,20 @@ def produce_inundated_branch_tif(
 
         # Use numpy.where operation to reclassify rem_path on the condition that the pixel values
         #   are <= to hand_stage and the catchments value is in the hydroid_list.
+
+        is_alaska = str(huc)[:2] == '19'
+
+        output_dtype = 'uint8' if is_alaska else 'int16'
+
         reclass_rem_array = np.where((rem_array <= hand_stage) & (rem_array != rem_src.nodata), 1, 0).astype(
-            'int16'
+            output_dtype
         )
 
         # The catchment_array has hydroid that have had the first 4 chars cut off
         # we need to the same for the hydroid's from the hydroid_list
         clipped_hydroid_list = []
         for i in hydroid_list:
-            clipped_str = str(i)[-4:]
+            clipped_str = str(i) if is_alaska else str(i)[-4:]
             clipped_hydroid_list.append(int(clipped_str))
 
         # Create a mask of the catchments_array where the values are in the clipped_hydroid_list
@@ -341,12 +348,12 @@ def produce_inundated_branch_tif(
         # Create a target array where the catchments_array is not nodata and the hydroid_mask is True
         target_catchments_array = np.where(
             ((hydroid_mask == True) & (catchments_array != catchments_src.nodata)), 1, 0
-        ).astype('int16')
+        ).astype(output_dtype)
 
         # Mask the reclass_rem_array with the target_catchments_array
         masked_reclass_rem_array = np.where(
             ((reclass_rem_array >= 1) & (target_catchments_array >= 1)), 1, 0
-        ).astype('int16')
+        ).astype(output_dtype)
 
         ## Debugging:
 

tools/convert_to_int16.py

@@ -1,14 +1,16 @@
 #!/usr/bin/env python3
-
 import argparse
 import os
+import sys
 import traceback
 from glob import glob
 
 import numpy as np
 import rioxarray as rxr
 import xarray as xr
 
+from utils.fim_enums import FIM_exit_codes
+
 
 def convert_to_int16(branch_dir: str):
     """
@@ -31,23 +33,17 @@ def convert_to_int16(branch_dir: str):
 
     # Iterate through each pair of gw catchments and rems
     for c, r in zip(catchments, rems):
-        rem = rxr.open_rasterio(r)
 
-        # Save original as another file to be deleted by deny list or saved
-        rem.rio.to_raster(r.replace('.tif', '_float32.tif'), compress="LZW", tiled=True)
-        nodata, crs = rem.rio.nodata, rem.rio.crs
-
-        # Preserve the second highest possible number for int16, use the highest number for nodata
-        rem.data = xr.where(rem > 32.766, 32.766, rem)
-        rem.data = xr.where(rem >= 0, np.round(rem * 1000), 32767)
-
-        rem = rem.astype(np.int16)
-        rem.rio.write_nodata(32767, inplace=True)
-        rem.rio.write_crs(crs, inplace=True)
+        catchment = rxr.open_rasterio(c)
 
-        rem.rio.to_raster(r, dtype=np.int16, compress="LZW", tiled=True)
+        # Check if converting data is possible
+        if (np.unique(catchment).shape[0] > 32766) | (len(str(int(np.max(catchment)))) > 8):
 
-        catchment = rxr.open_rasterio(c)
+            print(
+                "Catchment raster has either more than 32766 unique HydroIDs or has HydroIDs with more",
+                "than 8 digits.  Please adjust data accordingly before running Int16 data conversion.",
+            )
+            sys.exit(FIM_exit_codes.CANNOT_CONVERT_HYDROIDS_TO_INT16.value)
 
         if hydroid_prefix is None:
             hydroid_prefix = str(int(np.floor(catchment.max() / 10000)))
@@ -69,6 +65,22 @@ def convert_to_int16(branch_dir: str):
             with open(hydroid_prefix_path, 'w') as file:
                 file.write(hydroid_prefix)
 
+        rem = rxr.open_rasterio(r)
+
+        # Save original as another file to be deleted by deny list or saved
+        rem.rio.to_raster(r.replace('.tif', '_float32.tif'), compress="LZW", tiled=True)
+        nodata, crs = rem.rio.nodata, rem.rio.crs
+
+        # Preserve the second highest possible number for int16, use the highest number for nodata
+        rem.data = xr.where(rem > 32.766, 32.766, rem)
+        rem.data = xr.where(rem >= 0, np.round(rem * 1000), 32767)
+
+        rem = rem.astype(np.int16)
+        rem.rio.write_nodata(32767, inplace=True)
+        rem.rio.write_crs(crs, inplace=True)
+
+        rem.rio.to_raster(r, dtype=np.int16, compress="LZW", tiled=True)
+
 
 if __name__ == "__main__":