[1pt] Update stage-based CatFIM to prevent overflooding if stage + elevation is provided in lieu of stage

docs/CHANGELOG.md

@@ -1,6 +1,16 @@
 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.5.13.9 - 2025-01-24 - [PR#1399](https://github.com/NOAA-OWP/inundation-mapping/pull/1399)
+
+This update improves stage-based CatFIM by detecting and correcting instances where the stage value provided in the WRDS database is actually stage + elevation (which is actually water surface elevation and, uncaught, causes overflooding). 
+
+### Changes
+- `inundation-mapping/tools/catfim/generate_categorical_fim.py`: Added an update to detect and fix cases where WSE is provided in lieu of stage. Added `uncorrected_stage` and `is_interval` columns to output CSV.
+- `inundation-mapping/tools/catfim/generate_categorical_fim_mapping.py`: Added update to facilitate the new `is_interval` column.
+
+<br/><br/>
+
 ## v4.5.13.8 - 2025-01-24 - [PR#1405](https://github.com/NOAA-OWP/inundation-mapping/pull/1405)
 
 Removing the references to lid_to_run from CatFIM in order to keep the CatFIM scripts cleaner.  

tools/catfim/generate_categorical_fim.py

@@ -737,6 +737,29 @@ def iterate_through_huc_stage_based(
                 if not os.path.exists(mapping_lid_directory):
                     os.mkdir(mapping_lid_directory)
 
+                # Check whether stage value is actually a WSE value, and fix if needed:
+                # Get lowest stage value
+                lowest_stage_val = stage_values_df['stage_value'].min()
+
+                maximum_stage_threshold = 250  # TODO: Move to a variables file?
+
+                # Make an "rfc_stage" column for better documentation which shows the original
+                # uncorrect WRDS value before we adjsuted it for inundation
+                stage_values_df['rfc_stage'] = stage_values_df['stage_value']
+
+                # Stage value is larger than the elevation value AND greater than the
+                # maximum stage threshold, subtract the elev from the "stage" value
+                # to get the actual stage
+
+                if (lowest_stage_val > lid_altitude) and (lowest_stage_val > maximum_stage_threshold):
+                    stage_values_df['stage_value'] = stage_values_df['stage_value'] - lid_altitude
+                    MP_LOG.lprint(
+                        f"{huc_lid_id}: Lowest stage val > elev and higher than max stage thresh. Subtracted elev from stage vals to fix."
+                    )
+
+                # +++++++++++++++++++++++++++++
+                # This section is for inundating stages and intervals come later
+
                 # At this point we have at least one valid stage/category
                 # cyle through on the stages that are valid
                 # This are not interval values
@@ -763,7 +786,6 @@ def iterate_through_huc_stage_based(
                     # These are the up to 5 magnitudes being inundated at their stage value
                     (messages, hand_stage, datum_adj_wse, datum_adj_wse_m) = produce_stage_based_lid_tifs(
                         stage_value,
-                        False,
                         datum_adj_ft,
                         branch_dir,
                         lid_usgs_elev,
@@ -823,6 +845,9 @@ def iterate_through_huc_stage_based(
 
                 # MP_LOG.trace(f"non_rec_stage_values_df is {non_rec_stage_values_df}")
 
+                # +++++++++++++++++++++++++++++
+                # Creating interval tifs (if applicable)
+
                 # We already inundated and created files for the specific stages just not the intervals
                 # Make list of interval recs to be created
                 interval_list = []  # might stay empty
@@ -854,7 +879,6 @@ def iterate_through_huc_stage_based(
                                 executor.submit(
                                     produce_stage_based_lid_tifs,
                                     interval_stage_value,
-                                    True,
                                     datum_adj_ft,
                                     branch_dir,
                                     lid_usgs_elev,
@@ -912,7 +936,9 @@ def iterate_through_huc_stage_based(
                 # for threshold in categories:  (threshold and category are somewhat interchangeable)
                 # some may have failed inundation, which we will rectify later
                 MP_LOG.trace(f"{huc_lid_id}: updating threshhold values")
+
                 for threshold in valid_stage_names:
+
                     try:
 
                         # we don't know if the magnitude/stage can be mapped yes it hasn't been inundated
@@ -929,7 +955,12 @@ def iterate_through_huc_stage_based(
                                 'q': flows[threshold],
                                 'q_uni': flows['units'],
                                 'q_src': flows['source'],
-                                'stage': thresholds[threshold],
+                                'rfs_stage': stage_values_df.loc[stage_values_df['stage_name'] == threshold][
+                                    'rfc_stage'
+                                ],
+                                'stage': stage_values_df.loc[stage_values_df['stage_name'] == threshold][
+                                    'stage_value'
+                                ],
                                 'stage_uni': thresholds['units'],
                                 's_src': thresholds['source'],
                                 'wrds_time': thresholds['wrds_timestamp'],
@@ -1131,7 +1162,7 @@ def __calc_stage_intervals(non_rec_stage_values_df, past_major_interval_cap, huc
                 # MP_LOG.trace(f"{huc_lid_id}: Added interval value of {int_val}")
                 stage_values_claimed.append(int_val)
 
-    MP_LOG.lprint(f"{huc_lid_id} interval recs are {interval_recs}")
+    # MP_LOG.lprint(f"{huc_lid_id} interval recs are {interval_recs}")
 
     return interval_recs
 

tools/catfim/generate_categorical_fim_mapping.py

@@ -43,7 +43,6 @@
 # we will use an MP object either way
 def produce_stage_based_lid_tifs(
     stage_val,
-    is_interval_stage,
     datum_adj_ft,
     branch_dir,
     lid_usgs_elev,
@@ -665,10 +664,12 @@ def post_process_huc(
 
                     # careful. ft can be part of the site name, so only check part 3
                     interval_stage = None
+                    is_interval = False
                     if len(file_name_parts) >= 3 and "fti" in file_name_parts[2]:
                         try:
                             stage_val = file_name_parts[2].replace("fti", "")
                             interval_stage = float(stage_val)
+                            is_interval = True
                         except ValueError:
                             interval_stage = None
                             MP_LOG.error(
@@ -686,6 +687,7 @@ def post_process_huc(
                         magnitude,
                         nws_lid_attributes_filename,
                         interval_stage,
+                        is_interval,
                         parent_log_output_file,
                         child_log_file_prefix,
                     )
@@ -862,6 +864,7 @@ def reformat_inundation_maps(
     magnitude,
     nws_lid_attributes_filename,
     interval_stage,
+    is_interval,
     parent_log_output_file,
     child_log_file_prefix,
 ):
@@ -914,6 +917,7 @@ def reformat_inundation_maps(
         extent_poly_diss['version'] = fim_version
         extent_poly_diss['huc'] = huc
         extent_poly_diss['interval_stage'] = interval_stage
+        extent_poly_diss['is_interval'] = is_interval
 
         # Project to Web Mercator
         extent_poly_diss = extent_poly_diss.to_crs(VIZ_PROJECTION)
@@ -931,6 +935,9 @@ def reformat_inundation_maps(
         # already has an ahps_lid column which we want and not the nws_lid column
         extent_poly_diss = extent_poly_diss.drop(columns='nws_lid')
 
+        # Remove uncorrected stage from interval rows (to decrease potential for confusion)
+        extent_poly_diss.loc[extent_poly_diss['is_interval'] == True, 'stage_uncorrected'] = None
+
         # Save dissolved multipolygon
         handle = os.path.split(tif_to_process)[1].replace('.tif', '')
         diss_extent_filename = os.path.join(gpkg_dir, f"{huc}_{handle}.gpkg")