Model Validation with USGS Gauge Data

This example demonstrates how to validate HEC-RAS model results against observed USGS gauge data. We'll compare modeled flows from the Bald Eagle Creek model to observed data from a downstream gauge and calculate comprehensive validation metrics.

Workflow Overview

1. Run HEC-RAS model (or use existing results)
2. Extract modeled timeseries at location near gauge
3. Query observed USGS data from downstream gauge
4. Align timeseries for direct comparison
5. Calculate validation metrics (NSE, KGE, peak error, RMSE, bias)
6. Create comparison plots (timeseries, scatter, residuals)
7. Assess model performance and identify calibration needs

Use Case: Model Performance Assessment

This workflow enables: - Quantitative model performance evaluation - Identification of systematic biases - Guidance for calibration efforts - Reporting model accuracy for stakeholders

Setup and Imports

Python

# =============================================================================
# DEVELOPMENT MODE TOGGLE
# =============================================================================
USE_LOCAL_SOURCE = False  # <-- TOGGLE THIS

if USE_LOCAL_SOURCE:
    import sys
    from pathlib import Path
    local_path = str(Path.cwd().parent)
    if local_path not in sys.path:
        sys.path.insert(0, local_path)
    print(f"📁 LOCAL SOURCE MODE: Loading from {local_path}/ras_commander")
else:
    print("📦 PIP PACKAGE MODE: Loading installed ras-commander")

# Import ras-commander
from ras_commander import init_ras_project, RasExamples, ras, RasCmdr
from ras_commander.hdf import HdfResultsPlan
from ras_commander.usgs import (
    get_gauge_metadata,
    retrieve_flow_data,
    align_timeseries
)
from ras_commander.usgs.metrics import (
    nash_sutcliffe_efficiency,
    kling_gupta_efficiency,
    calculate_peak_error,
    calculate_all_metrics
)
from ras_commander.usgs.visualization import (
    plot_timeseries_comparison,
    plot_scatter_comparison,
    plot_residuals
)

# Additional imports
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime, timedelta

# Verify which version loaded
import ras_commander
print(f"✓ Loaded: {ras_commander.__file__}")
print("✓ Imports successful")

Model Validation Metrics

After calculating calibration/validation metrics:

Target Metrics (Moriasi et al. 2007):

Metric	Excellent	Good	Satisfactory	Unsatisfactory
NSE	>0.75	0.65-0.75	0.50-0.65	<0.50
PBIAS	<±10%	±10-15%	±15-25%	>±25%
R²	>0.85	0.75-0.85	0.60-0.75	<0.60

Verification Checklist: - [ ] NSE >0.50 (satisfactory minimum) - [ ] PBIAS <±25% (volume error acceptable) - [ ] Visual inspection confirms hydrograph shape match - [ ] Peak flow timing within 1-2 hours - [ ] Recession limb shape matches observed

Documentation: - Plot observed vs simulated hydrographs (include in submittal) - Table of metrics by validation event - Narrative explaining any poor-fitting events (ice, debris, etc.)

References: - Moriasi et al. (2007) Model Evaluation Guidelines - Moriasi et al. (2015) Updated Guidelines - ASCE/EWRI (2017) Flood Forecasting & Warning Systems

Parameters

Configure these values to customize the notebook for your project.

Python

# =============================================================================
# PARAMETERS - Edit these to customize the notebook
# =============================================================================
from pathlib import Path

# Project Configuration
PROJECT_NAME = "Muncie"           # Example project to extract
RAS_VERSION = "7.0"               # HEC-RAS version (6.3, 6.5, 6.6, etc.)
SUFFIX = "914"                    # Notebook identifier for extraction

# USGS Configuration
USGS_SITE = "01547200"            # USGS gauge site number
START_DATE = "2020-01-01"         # Data start date
END_DATE = "2020-12-31"           # Data end date
ONLINE = True                     # Enable network requests

print(f"Configured for project: {PROJECT_NAME}")
print(f"RAS Version: {RAS_VERSION}")
print(f"Extraction suffix: {SUFFIX}")

1. Initialize Project and Check Results

Python

# Extract and initialize Bald Eagle Creek project
print("Extracting Bald Eagle Creek project...")
project_paths = RasExamples.extract_project(["Balde Eagle Creek", "BaldEagleCrkMulti2D"], suffix=SUFFIX)
project_path = project_paths[0]  # Use 1D model
print(f"Using: {project_path}\n")

init_ras_project(project_path, RAS_VERSION)
print(f"Project initialized: {ras.project_folder}")
print(f"Plans: {len(ras.plan_df)}")

Python

# Check if results exist
plan_df = ras.plan_df
plan_01 = plan_df[plan_df['plan_number'] == '01'].iloc[0]

results_path = plan_01.get('results_path') or plan_01.get('hdf_path')
has_results = Path(results_path).exists() if results_path else False

print(f"\nPlan 01 Results:")
print(f"  Results file: {Path(results_path).name if results_path else 'Not found'}")
print(f"  Exists: {has_results}")

if not has_results:
    print(f"\n⚠ No results found - will run model first")
else:
    print(f"\n✓ Results exist - can proceed with validation")

2. Run Model (if needed)

Python

# Run model if results don't exist
if not has_results:
    print("Running HEC-RAS Plan 01...")
    print("This may take several minutes...\n")

    try:
        RasCmdr.compute_plan(
            plan_number='01',
            ras_object=ras,
            num_cores=8
        )
        print("\n✓ Model execution complete")

        # Reinitialize to load new results
        init_ras_project(project_path, RAS_VERSION)
        plan_01 = ras.plan_df[ras.plan_df['plan_number'] == '01'].iloc[0]
        results_path = plan_01.get('results_path') or plan_01.get('hdf_path')

    except Exception as e:
        print(f"\n✗ Error running model: {e}")
        print("This example requires HEC-RAS to be installed.")
        results_path = None
else:
    print("Using existing model results")

3. Extract Modeled Timeseries

We'll extract modeled flow at a cross section near the downstream gauge location.

Validation Gauge: USGS-01548005 (Bald Eagle Creek near Beech Creek Station) - Location: ~10 miles downstream of Lock Haven - Drainage: 562 sq mi - Model XS: We'll use a cross section in the downstream reach

Python

# Open results HDF
if results_path and Path(results_path).exists():
    print(f"Opening results: {Path(results_path).name}")
    hdf = HdfResultsPlan(results_path)

    # Get available cross sections
    xs_list = hdf.get_cross_section_list()
    print(f"\nAvailable cross sections: {len(xs_list)}")
    print(f"River stations range: {xs_list[0]} to {xs_list[-1]}")

    # Select a downstream cross section near the gauge
    # USGS-01548005 is ~10 miles downstream, likely around river station 80000-90000
    target_station = 81500  # Near the gate structure, downstream location

    # Find closest cross section
    closest_xs = min(xs_list, key=lambda x: abs(float(x) - target_station))

    print(f"\nValidation location:")
    print(f"  Target station: {target_station}")
    print(f"  Closest XS: {closest_xs}")
    print(f"  This is downstream of Lock Haven, near USGS-01548005")

else:
    print("Cannot extract modeled data - no results available")
    closest_xs = None

Python

# Extract modeled flow timeseries
if closest_xs:
    print(f"Extracting modeled flow at XS {closest_xs}...")

    # Get flow timeseries
    flow_ts = hdf.get_timeseries(
        river='Bald Eagle',
        reach='Loc Hav',
        rs=closest_xs,
        variable='Flow'
    )

    # Convert to DataFrame
    modeled_df = pd.DataFrame({
        'datetime': flow_ts['datetime'],
        'value': flow_ts['values']
    })

    print(f"\nModeled flow statistics:")
    print(f"  Records: {len(modeled_df)}")
    print(f"  Time range: {modeled_df['datetime'].min()} to {modeled_df['datetime'].max()}")
    print(f"  Mean flow: {modeled_df['value'].mean():.1f} cfs")
    print(f"  Peak flow: {modeled_df['value'].max():.1f} cfs")
    print(f"  Min flow: {modeled_df['value'].min():.1f} cfs")
else:
    print("Cannot extract modeled timeseries - no results")
    modeled_df = None

4. Query Observed USGS Data

Retrieve observed flow data from the downstream gauge for the same time period.

Python

# Gauge configuration
downstream_gauge = "01548005"  # Bald Eagle Creek near Beech Creek Station

print(f"Querying USGS-{downstream_gauge}...")
try:
    metadata = get_gauge_metadata(downstream_gauge)
    print(f"\nGauge: USGS-{downstream_gauge}")
    print(f"Name: {metadata['station_name']}")
    print(f"Location: ({metadata['latitude']:.4f}, {metadata['longitude']:.4f})")
    print(f"Drainage Area: {metadata['drainage_area_sqmi']} sq mi")
except Exception as e:
    print(f"Note: Metadata query had issues: {e}")
    print(f"Using known values: USGS-{downstream_gauge}, Beech Creek Station, PA, 562 sq mi")

Python

# Get simulation period from model
if modeled_df is not None:
    start_date = modeled_df['datetime'].min()
    end_date = modeled_df['datetime'].max()

    print(f"\nQuerying observed data for model period:")
    print(f"  Start: {start_date}")
    print(f"  End: {end_date}")

    # Retrieve observed data
    observed_df = retrieve_flow_data(
        site_id=downstream_gauge,
        start_date=start_date.strftime('%Y-%m-%d'),
        end_date=end_date.strftime('%Y-%m-%d'),
        interval='instantaneous'  # Hourly data
    )

    print(f"\nObserved flow statistics:")
    print(f"  Records: {len(observed_df)}")
    print(f"  Time range: {observed_df['datetime'].min()} to {observed_df['datetime'].max()}")
    print(f"  Mean flow: {observed_df['value'].mean():.1f} cfs")
    print(f"  Peak flow: {observed_df['value'].max():.1f} cfs")
    print(f"  Min flow: {observed_df['value'].min():.1f} cfs")

else:
    print("Cannot query observed data - no model results")
    observed_df = None

5. Align Timeseries for Comparison

Synchronize modeled and observed data to matching timestamps.

Python

# Align timeseries
if modeled_df is not None and observed_df is not None:
    print("Aligning modeled and observed timeseries...")

    aligned_df = align_timeseries(
        modeled=modeled_df,
        observed=observed_df,
        tolerance_minutes=30  # Allow 30-minute tolerance for timestamp matching
    )

    print(f"\nAligned timeseries:")
    print(f"  Matched records: {len(aligned_df)}")
    print(f"  Time range: {aligned_df['datetime'].min()} to {aligned_df['datetime'].max()}")
    print(f"  Data completeness: {len(aligned_df) / len(modeled_df) * 100:.1f}%")

    print(f"\nAligned statistics:")
    print(f"  Modeled mean: {aligned_df['modeled'].mean():.1f} cfs")
    print(f"  Observed mean: {aligned_df['observed'].mean():.1f} cfs")
    print(f"  Difference: {aligned_df['modeled'].mean() - aligned_df['observed'].mean():.1f} cfs")

else:
    print("Cannot align timeseries - missing data")
    aligned_df = None

6. Calculate Validation Metrics

Compute comprehensive performance metrics for model assessment.

Python

# Calculate all validation metrics
if aligned_df is not None and len(aligned_df) > 0:
    print("Calculating validation metrics...\n")

    metrics = calculate_all_metrics(
        observed=aligned_df['observed'].values,
        modeled=aligned_df['modeled'].values
    )

    print("=" * 70)
    print("MODEL PERFORMANCE METRICS")
    print("=" * 70)

    print(f"\nOverall Performance:")
    print(f"  Nash-Sutcliffe Efficiency (NSE): {metrics['nse']:.3f}")
    print(f"    Interpretation: {_interpret_nse(metrics['nse'])}")
    print(f"  Kling-Gupta Efficiency (KGE): {metrics['kge']:.3f}")
    print(f"    Components: r={metrics['kge_components']['correlation']:.3f}, "
          f"β={metrics['kge_components']['bias_ratio']:.3f}, "
          f"γ={metrics['kge_components']['variability_ratio']:.3f}")

    print(f"\nError Metrics:")
    print(f"  Root Mean Square Error (RMSE): {metrics['rmse']:.1f} cfs")
    print(f"  Mean Absolute Error (MAE): {metrics['mae']:.1f} cfs")
    print(f"  Mean Bias: {metrics['bias']:.1f} cfs")
    print(f"    Model {'overestimates' if metrics['bias'] > 0 else 'underestimates'} by {abs(metrics['bias']):.1f} cfs on average")

    print(f"\nPeak Flow Analysis:")
    peak_metrics = calculate_peak_error(
        observed=aligned_df['observed'].values,
        modeled=aligned_df['modeled'].values,
        timestamps=aligned_df['datetime'].values
    )
    print(f"  Observed peak: {peak_metrics['observed_peak']:.1f} cfs at {peak_metrics['observed_peak_time']}")
    print(f"  Modeled peak: {peak_metrics['modeled_peak']:.1f} cfs at {peak_metrics['modeled_peak_time']}")
    print(f"  Peak error: {peak_metrics['peak_error_percent']:.1f}%")
    print(f"  Timing error: {peak_metrics['timing_error_hours']:.1f} hours")

    print(f"\nCorrelation:")
    print(f"  Pearson R²: {metrics['r_squared']:.3f}")
    print(f"  Percent Bias (PBIAS): {metrics['pbias']:.1f}%")

    print("=" * 70)

else:
    print("Cannot calculate metrics - no aligned data")
    metrics = None

# Helper function for NSE interpretation
def _interpret_nse(nse):
    if nse > 0.75:
        return "Excellent (> 0.75)"
    elif nse > 0.65:
        return "Very Good (0.65-0.75)"
    elif nse > 0.50:
        return "Good (0.50-0.65)"
    elif nse > 0.40:
        return "Satisfactory (0.40-0.50)"
    else:
        return "Unsatisfactory (< 0.40)"

7. Create Comparison Plots

Generate publication-quality visualizations of model performance.

Python

# Timeseries comparison plot
if aligned_df is not None and metrics is not None:
    print("Creating timeseries comparison plot...")

    plot_timeseries_comparison(
        observed=aligned_df['observed'].values,
        modeled=aligned_df['modeled'].values,
        timestamps=aligned_df['datetime'].values,
        title=f"Model Validation: USGS-{downstream_gauge} vs HEC-RAS",
        observed_label="USGS Observed",
        modeled_label="HEC-RAS Modeled",
        show_metrics=True,
        metrics_dict=metrics
    )

    plt.tight_layout()
    plt.show()

Python

# Scatter plot with 1:1 line
if aligned_df is not None:
    print("Creating scatter comparison plot...")

    plot_scatter_comparison(
        observed=aligned_df['observed'].values,
        modeled=aligned_df['modeled'].values,
        title="Modeled vs Observed Flow",
        xlabel="Observed Flow (cfs)",
        ylabel="Modeled Flow (cfs)",
        show_metrics=True,
        metrics_dict=metrics
    )

    plt.tight_layout()
    plt.show()

Python

# Residual analysis (4-panel diagnostic plot)
if aligned_df is not None:
    print("Creating residual diagnostic plots...")

    plot_residuals(
        observed=aligned_df['observed'].values,
        modeled=aligned_df['modeled'].values,
        timestamps=aligned_df['datetime'].values
    )

    plt.tight_layout()
    plt.show()

8. Performance Assessment and Recommendations

Python

# Generate performance assessment
if metrics is not None:
    print("=" * 70)
    print("PERFORMANCE ASSESSMENT")
    print("=" * 70)

    nse = metrics['nse']
    kge = metrics['kge']
    pbias = metrics['pbias']

    print(f"\nOverall Model Rating:")
    if nse > 0.75 and abs(pbias) < 10:
        rating = "EXCELLENT"
        color = "green"
    elif nse > 0.65 and abs(pbias) < 15:
        rating = "VERY GOOD"
        color = "lightgreen"
    elif nse > 0.50 and abs(pbias) < 25:
        rating = "GOOD"
        color = "yellow"
    elif nse > 0.40:
        rating = "SATISFACTORY"
        color = "orange"
    else:
        rating = "UNSATISFACTORY"
        color = "red"

    print(f"  {rating}")
    print(f"  NSE: {nse:.3f}, KGE: {kge:.3f}, PBIAS: {pbias:.1f}%")

    print(f"\nKey Findings:")

    # Bias assessment
    if abs(pbias) < 10:
        print(f"  ✓ Low bias ({pbias:.1f}%) - Model captures flow magnitude well")
    elif abs(pbias) < 25:
        print(f"  ⚠ Moderate bias ({pbias:.1f}%) - Consider adjusting Manning's n or BC scaling")
    else:
        print(f"  ✗ High bias ({pbias:.1f}%) - Significant systematic error present")

    # Timing assessment
    if 'timing_error_hours' in peak_metrics:
        timing_error = abs(peak_metrics['timing_error_hours'])
        if timing_error < 3:
            print(f"  ✓ Good timing ({timing_error:.1f}h error) - Peak timing well captured")
        elif timing_error < 6:
            print(f"  ⚠ Moderate timing error ({timing_error:.1f}h) - Check wave travel time")
        else:
            print(f"  ✗ Poor timing ({timing_error:.1f}h) - Review time step and boundary timing")

    # Peak magnitude assessment
    if 'peak_error_percent' in peak_metrics:
        peak_error = abs(peak_metrics['peak_error_percent'])
        if peak_error < 15:
            print(f"  ✓ Good peak match ({peak_error:.1f}% error) - Peak flows accurate")
        elif peak_error < 25:
            print(f"  ⚠ Moderate peak error ({peak_error:.1f}%) - Review roughness calibration")
        else:
            print(f"  ✗ Large peak error ({peak_error:.1f}%) - Significant peak flow mismatch")

    print(f"\nCalibration Recommendations:")

    if pbias > 15:
        print(f"  1. Model overestimates - Increase Manning's n or reduce BC flows")
    elif pbias < -15:
        print(f"  1. Model underestimates - Decrease Manning's n or increase BC flows")
    else:
        print(f"  1. Flow magnitude is well calibrated")

    if nse < 0.65:
        print(f"  2. Improve NSE by reducing systematic errors (check geometry, roughness)")

    if metrics['kge_components']['correlation'] < 0.85:
        print(f"  3. Low correlation - Check boundary condition timing and shape")

    print("\n" + "=" * 70)

Summary: Model Validation Workflow

This example demonstrated the complete workflow for validating HEC-RAS results against USGS gauge observations:

Workflow Steps Completed

1. ✓ Ran HEC-RAS model (or used existing results)
2. ✓ Extracted modeled flow at cross section near downstream gauge
3. ✓ Queried USGS observations from gauge 01548005 (Beech Creek Station)
4. ✓ Aligned timeseries for direct comparison
5. ✓ Calculated validation metrics (NSE, KGE, RMSE, peak error, etc.)
6. ✓ Created comparison plots (timeseries, scatter, residuals)
7. ✓ Assessed performance and provided calibration guidance

Key Functions Used

• retrieve_flow_data() - Query USGS observations
• align_timeseries() - Synchronize modeled and observed data
• calculate_all_metrics() - Comprehensive validation suite
• nash_sutcliffe_efficiency() - NSE metric
• kling_gupta_efficiency() - KGE metric with components
• calculate_peak_error() - Peak flow and timing analysis
• plot_timeseries_comparison() - Publication-quality timeseries plots
• plot_scatter_comparison() - 1:1 scatter plots with statistics
• plot_residuals() - 4-panel diagnostic plots

Validation Metrics Guide

Nash-Sutcliffe Efficiency (NSE): - > 0.75: Excellent - 0.65-0.75: Very Good - 0.50-0.65: Good - 0.40-0.50: Satisfactory - < 0.40: Unsatisfactory

Percent Bias (PBIAS): - < ±10%: Very Good - ±10-15%: Good - ±15-25%: Satisfactory - > ±25%: Unsatisfactory

Kling-Gupta Efficiency (KGE): - Balances correlation, bias, and variability - Values near 1.0 indicate excellent performance - Components help identify specific deficiencies

Applications

This workflow enables: - Quantitative model performance reporting - Systematic calibration guidance - Identification of bias and timing errors - Stakeholder communication of model accuracy

Related Examples

• Example 29: Complete USGS integration workflow
• 912_usgs_real_time_monitoring.ipynb: Real-time monitoring functions
• 913_bc_generation_from_live_gauge.ipynb: BC generation from live gauge data