[5pt] PR: Bug fixes for Sierra test

docs/CHANGELOG.md

@@ -1,6 +1,13 @@
 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.x.x - 2025-01-02 - [PR#1388]https://github.com/NOAA-OWP/inundation-mapping/pull/1388)
+
+Fixed Sierra test bugs to draw the vertical lines.
+
+### Changes
+
+- `tools/rating_curve_comparison.py`
 
 
 ## v4.5.13.4 - 2024-01-03 - [PR#1382](https://github.com/NOAA-OWP/inundation-mapping/pull/1382)

tools/rating_curve_comparison.py

@@ -255,7 +255,9 @@ def generate_rating_curve_metrics(args):
                 str_order = np.unique(usgs_rc.order_).item()
                 feature_id = str(gage.feature_id)
 
-                usgs_pred_elev = get_reccur_intervals(usgs_rc, usgs_crosswalk, nwm_recurr_intervals_all)
+                usgs_pred_elev, feature_index = get_reccur_intervals(
+                    usgs_rc, usgs_crosswalk, nwm_recurr_intervals_all
+                )
 
                 # Handle sites missing data
                 if len(usgs_pred_elev) < 1:
@@ -282,7 +284,10 @@ def generate_rating_curve_metrics(args):
                     )
                     continue
 
-                fim_pred_elev = get_reccur_intervals(fim_rc, usgs_crosswalk, nwm_recurr_intervals_all)
+                if feature_index is not None:
+                    fim_pred_elev = get_reccur_intervals_fim(
+                        fim_rc, usgs_crosswalk, nwm_recurr_intervals_all, feature_index
+                    )
 
                 # Handle sites missing data
                 if len(fim_pred_elev) < 1:
@@ -942,14 +947,14 @@ def generate_rc_and_rem_plots(rc, plot_filename, recurr_data_table, branches_fol
     plt.close()
 
 
-def get_reccur_intervals(site_rc, usgs_crosswalk, nwm_recurr_intervals):
+def get_reccur_intervals_fim(site_rc, usgs_crosswalk, nwm_recurr_intervals, feature_index):
     usgs_site = site_rc.merge(usgs_crosswalk, on="location_id")
     nwm_ids = len(usgs_site.feature_id.drop_duplicates())
 
     if nwm_ids > 0:
         try:
             nwm_recurr_intervals = nwm_recurr_intervals.copy().loc[
-                nwm_recurr_intervals.feature_id == usgs_site.feature_id.drop_duplicates().item()
+                nwm_recurr_intervals.feature_id == usgs_site.feature_id.drop_duplicates().loc[feature_index]
             ]
             nwm_recurr_intervals['pred_elev'] = np.interp(
                 nwm_recurr_intervals.discharge_cfs.values,
@@ -972,6 +977,55 @@ def get_reccur_intervals(site_rc, usgs_crosswalk, nwm_recurr_intervals):
         return []
 
 
+def get_reccur_intervals(site_rc, usgs_crosswalk, nwm_recurr_intervals):
+    usgs_site = site_rc.merge(usgs_crosswalk, on="location_id")
+    nwm_ids = len(usgs_site.feature_id.drop_duplicates())
+    min_discharge = site_rc.loc[(site_rc.source == 'USGS')].discharge_cfs.min()
+    max_discharge = site_rc.loc[(site_rc.source == 'USGS')].discharge_cfs.max()
+    discharge_range = max_discharge - min_discharge
+
+    if nwm_ids > 0:
+        try:
+            filtered = nwm_recurr_intervals.copy().loc[
+                nwm_recurr_intervals.feature_id == usgs_site.feature_id.drop_duplicates().iloc[0]
+            ]
+            min_q_recurr = filtered.discharge_cfs.min()
+            max_q_recurr = filtered.discharge_cfs.max()
+            spread_q = max_q_recurr - min_q_recurr
+            ratio = spread_q / discharge_range
+            # If there is only one feature_id for each location_id
+            if nwm_ids == 1:
+                nwm_recurr_intervals = filtered
+                feature_index = 0
+
+            # If there is more one feature_id for each location_id
+            else:
+                if ratio > 0.1:
+                    nwm_recurr_intervals = filtered
+                    feature_index = 0
+                else:
+                    nwm_recurr_intervals = nwm_recurr_intervals.copy().loc[
+                        nwm_recurr_intervals.feature_id == usgs_site.feature_id.drop_duplicates().iloc[1]
+                    ]
+                    feature_index = 1
+            nwm_recurr_intervals['pred_elev'] = np.interp(
+                nwm_recurr_intervals.discharge_cfs.values,
+                usgs_site['discharge_cfs'],
+                usgs_site['elevation_ft'],
+                left=np.nan,
+                right=np.nan,
+            )
+
+            return nwm_recurr_intervals, feature_index
+        except Exception as ex:
+            summary = traceback.StackSummary.extract(traceback.walk_stack(None))
+            print(summary, repr(ex))
+            return []
+
+    else:
+        return []
+
+
 def calculate_rc_stats_elev(rc, stat_groups=None):
     usgs_elev = "USGS"
     src_elev = "FIM"
@@ -1299,6 +1353,7 @@ def evaluate_results(sierra_results=[], labels=[], save_location=''):
     )  # using WARNING level to avoid benign? info messages ("Failed to auto identify EPSG: 7")
     format = '  %(message)s'
     log_dt_string = start_time.strftime("%Y_%m_%d-%H_%M_%S")
+    os.makedirs(output_dir, exist_ok=True)
     handlers = [
         logging.FileHandler(os.path.join(output_dir, f'rating_curve_comparison_{log_dt_string}.log')),
         logging.StreamHandler(),