Skip to content

Historical Event Validation with AORC Precipitation and USGS Gauges

This notebook demonstrates a comprehensive historical flood event validation workflow using: - AORC gridded precipitation (rain-on-grid on 2D mesh) - USGS gauge data for boundary conditions - Multiple validation points - HUC12 watershed coverage analysis

Event: December 24-25, 2020 storm (2.72 inches, largest 2020 event) Model: Bald Eagle Creek Multi-2D
Template: Plan 06 (gridded precipitation enabled)

Workflow Overview

  1. Extract and initialize HEC-RAS project (BaldEagleCrkMulti2D)
  2. Get project bounds and HUC12 watershed coverage analysis
  3. Download AORC precipitation data for storm event
  4. Retrieve USGS gauge data for boundary conditions
  5. Clone plan and configure for storm simulation
  6. Run HEC-RAS model
  7. Extract modeled results and compare with observed USGS data
  8. Calculate validation metrics and generate comparison plots
Python
# =============================================================================
# IMPORTS AND SETUP
# =============================================================================
from pathlib import Path
import sys

# Flexible imports for development vs installed package
try:
    from ras_commander import RasExamples, init_ras_project, RasCmdr, RasPlan, RasUnsteady, RasUtils, ras
    from ras_commander.hdf import HdfProject, HdfMesh, HdfResultsXsec, HdfResultsMesh
    from ras_commander.precip import PrecipAorc
    from ras_commander.usgs import (
        get_gauge_metadata,
        retrieve_flow_data,
        retrieve_stage_data,
        align_timeseries,
        calculate_all_metrics,
        plot_timeseries_comparison,
        plot_scatter_comparison,
        configure_rate_limit
    )
except ImportError:
    current_file = Path.cwd()
    parent_directory = current_file.parent
    sys.path.insert(0, str(parent_directory))
    from ras_commander import RasExamples, init_ras_project, RasCmdr, RasPlan, RasUnsteady, RasUtils, ras
    from ras_commander.hdf import HdfProject, HdfMesh, HdfResultsXsec, HdfResultsMesh
    from ras_commander.precip import PrecipAorc
    from ras_commander.usgs import (
        get_gauge_metadata,
        retrieve_flow_data,
        retrieve_stage_data,
        align_timeseries,
        calculate_all_metrics,
        plot_timeseries_comparison,
        plot_scatter_comparison,
        configure_rate_limit
    )

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.geometry import Point

print("Imports successful")
Text Only
Imports successful
Python
# =============================================================================
# API KEY SETUP
# =============================================================================
# USGS API key provides higher rate limits (5 req/sec vs 0.2 req/sec)
api_key_file = Path("usgs_api_key.txt")

if api_key_file.exists():
    usgs_api_key = api_key_file.read_text().strip()
    configure_rate_limit(requests_per_second=5.0)
    print(f"USGS API key loaded - configured rate limit: 5.0 req/sec")
else:
    usgs_api_key = None
    print("No API key file found - using public access (0.2 req/sec)")
    print("To get higher rate limits, create usgs_api_key.txt with your USGS API key")
Text Only
No API key file found - using public access (0.2 req/sec)
To get higher rate limits, create usgs_api_key.txt with your USGS API key
Python
# =============================================================================
# PARAMETERS
# =============================================================================
PROJECT_NAME = "BaldEagleCrkMulti2D"
TEMPLATE_PLAN = "06"  # Plan with gridded precipitation enabled
STORM_DATE = "20201224"  # December 24-25, 2020 storm

# Simulation period (48h warmup + event + 48h recession)
SIM_START = "2020-12-22"
SIM_END = "2020-12-27"

# USGS Gauges in Bald Eagle Creek watershed
# 01547500 - Bald Eagle Creek at Blanchard, PA (upstream BC)
# 01548005 - Beech Creek Station (validation point + downstream BC)
UPSTREAM_GAUGE = "01547500"
VALIDATION_GAUGE = "01548005"

# HEC-RAS version
RAS_VERSION = "7.0"

print(f"Project: {PROJECT_NAME}")
print(f"Template Plan: {TEMPLATE_PLAN}")
print(f"Storm Date: {STORM_DATE}")
print(f"Simulation Period: {SIM_START} to {SIM_END}")
print(f"Upstream BC Gauge: {UPSTREAM_GAUGE}")
print(f"Validation Gauge: {VALIDATION_GAUGE}")
Text Only
Project: BaldEagleCrkMulti2D
Template Plan: 06
Storm Date: 20201224
Simulation Period: 2020-12-22 to 2020-12-27
Upstream BC Gauge: 01547500
Validation Gauge: 01548005
Python
# =============================================================================
# EXTRACT AND INITIALIZE PROJECT
# =============================================================================
suffix = "914_historical"
expected_folder = Path.cwd() / "example_projects" / f"{PROJECT_NAME}_{suffix}"

# Recreate the project copy each run so cloned plans, HDFs, and precipitation
# metadata always reflect the current notebook cells rather than stale artifacts.
if expected_folder.exists():
    print(f"Removing existing project folder: {expected_folder}")
    if not RasUtils.remove_with_retry(expected_folder, ras_object=None):
        raise PermissionError(f"Unable to remove stale project folder: {expected_folder}")

project_folder = RasExamples.extract_project(PROJECT_NAME, suffix=suffix)
print(f"Project extracted to: {project_folder}")

# Initialize project
ras = init_ras_project(project_folder, RAS_VERSION)
print("Project initialized")

# Display available plans
print("\nAvailable plans:")
print(ras.plan_df[['plan_number', 'Plan Title']].to_string())
Text Only
Removing existing project folder: C:\Users\bill\.config\superpowers\worktrees\ras-commander-pr-docs\codex-pr251-timeseries-alignment\examples\example_projects\BaldEagleCrkMulti2D_914_historical


2026-07-08 21:38:37 - ras_commander.RasExamples - INFO - Successfully extracted project 'BaldEagleCrkMulti2D' to example_projects\BaldEagleCrkMulti2D_914_historical


Project extracted to: C:\Users\bill\.config\superpowers\worktrees\ras-commander-pr-docs\codex-pr251-timeseries-alignment\examples\example_projects\BaldEagleCrkMulti2D_914_historical


2026-07-08 21:38:38 - ras_commander.RasPrj - INFO - ras-commander v0.98.2 | An open-source project of CLB Engineering Corporation (https://clbengineering.com/) | Docs: https://rascommander.info | GitHub: https://github.com/gpt-cmdr/ras-commander


2026-07-08 21:38:38 - ras_commander.RasPrj - INFO - Project initialized: BaldEagleDamBrk


2026-07-08 21:38:38 - ras_commander.RasPrj - INFO - Using HEC-RAS version 7.0


2026-07-08 21:38:38 - ras_commander.RasPrj - INFO - 
═══════════════════════════════════════════════════════════════════════
ras-commander | HEC-RAS Automation Library
Docs: https://rascommander.info/
Repo: https://github.com/gpt-cmdr/ras-commander
LLM agents: https://rascommander.info/llms.txt
═══════════════════════════════════════════════════════════════════════

PROJECT DATAFRAMES (single source of truth — use these, not file globbing):
  ras.plan_df        Plans, HDF paths, geometry/flow associations
  ras.geom_df        Geometry files and HDF preprocessor paths
  ras.flow_df        Steady flow files
  ras.unsteady_df    Unsteady flow files and configurations
  ras.boundaries_df  Boundary conditions (type, name, location)
  ras.results_df     Lightweight HDF results summaries
  ras.rasmap_df      RASMapper layers, terrain, land cover paths

KEY APIS (static classes — call directly, never instantiate):
  Execution:    RasCmdr.compute_plan() / compute_parallel() / compute_test_mode()
  Plan Files:   RasPlan.clone_plan() / clone_geom() / set_geom()
  Unsteady:     RasUnsteady — IC/BC management, gate openings, precipitation
  Geometry:     GeomCrossSection, GeomBridge, GeomStorage, GeomLateral, GeomMesh
  HDF Results:  HdfResultsPlan.get_wse() / get_compute_messages()
                HdfResultsMesh.get_mesh_max_ws() / get_mesh_cells_timeseries()
                HdfMesh.get_mesh_cell_points()
  QA/QC:        RasCheck.run_check() / RasFixit (geometry repair)
  DSS:          RasDss.get_timeseries() / check_pathname()
  USGS:         UsgsGaugeSpatial, GaugeMatcher, RasUsgsBoundaryGeneration
  Precipitation: StormGenerator, Atlas14Storm, PrecipAorc, Atlas14Variance
  Terrain:      RasTerrain.create_terrain_hdf() / RasTerrainMod

MULTI-PROJECT: Pass ras_object= to all API calls when using local RasPrj instances.

EXAMPLES: 100+ notebooks in examples/ (100s=execution, 200s=geometry, 300s=unsteady,
  400s=HDF results, 500s=remote, 800s=QA/QC, 900s=data integration).
  Review relevant notebooks before assembling new workflows.

PLATFORM: Most HEC-RAS operations require Windows. Linux/Wine support for
  headless execution, data access, geometry modification, and preprocessing
  is available via RasProcess (HEC-RAS 6.6+). See ras_commander/RasProcess.py.
  Remote distributed execution: ras_commander/remote/ (PsExec, Docker, SSH, cloud).
═══════════════════════════════════════════════════════════════════════


Project initialized

Available plans:
   plan_number                               Plan Title
0           13                  PMF with Multi 2D Areas
1           15              1d-2D Dambreak Refined Grid
2           17                          2D to 1D No Dam
3           18                             2D to 2D Run
4           19                   SA to 2D Dam Break Run
5           03  Single 2D Area - Internal Dam Structure
6           04  SA to 2D Area Conn - 2D Levee Structure
7           02                 SA to Detailed 2D Breach
8           01             SA to Detailed 2D Breach FEQ
9           05          Single 2D area with Bridges FEQ
10          06            Gridded Precip - Infiltration
Python
# =============================================================================
# GET PROJECT BOUNDS FOR AORC DOWNLOAD
# =============================================================================
# Find geometry HDF file
geom_hdf_files = list(project_folder.glob("*.g*.hdf"))
print(f"Found geometry HDF files: {[f.name for f in geom_hdf_files]}")

# Use the most appropriate geometry file (typically highest number with 2D areas)
geom_hdf = None
for f in sorted(geom_hdf_files, reverse=True):
    if f.suffix == '.hdf' and '.g' in f.name:
        geom_hdf = f
        break

if geom_hdf is None:
    raise FileNotFoundError("No geometry HDF file found")

print(f"\nUsing geometry HDF: {geom_hdf.name}")

# Get project bounds in WGS84 for AORC download
# Use 50% buffer to capture upstream precipitation areas
bounds = HdfProject.get_project_bounds_latlon(
    geom_hdf, 
    buffer_percent=50.0,
    project_crs="EPSG:2271"  # NAD83 StatePlane Pennsylvania North, US feet
)

print(f"\nProject bounds (WGS84 with 50% buffer):")
print(f"  West:  {bounds[0]:.4f}")
print(f"  South: {bounds[1]:.4f}")
print(f"  East:  {bounds[2]:.4f}")
print(f"  North: {bounds[3]:.4f}")
Text Only
Found geometry HDF files: ['BaldEagleDamBrk.g01.hdf', 'BaldEagleDamBrk.g02.hdf', 'BaldEagleDamBrk.g03.hdf', 'BaldEagleDamBrk.g06.hdf', 'BaldEagleDamBrk.g08.hdf', 'BaldEagleDamBrk.g09.hdf', 'BaldEagleDamBrk.g10.hdf', 'BaldEagleDamBrk.g11.hdf', 'BaldEagleDamBrk.g12.hdf', 'BaldEagleDamBrk.g13.hdf']

Using geometry HDF: BaldEagleDamBrk.g13.hdf

Project bounds (WGS84 with 50% buffer):
  West:  -77.7085
  South: 41.0102
  East:  -77.2511
  North: 41.2197
Python
# =============================================================================
# DOWNLOAD HUC12 WATERSHED FOR COVERAGE ANALYSIS
# =============================================================================
try:
    from pygeohydro import WBD

    # Get mesh areas for centroid calculation
    mesh_areas = HdfMesh.get_mesh_areas(geom_hdf)
    print(f"Found {len(mesh_areas)} 2D flow areas")

    # Calculate centroid of 2D mesh in WGS84
    mesh_wgs84 = mesh_areas.to_crs("EPSG:4326")
    centroid = mesh_wgs84.geometry.unary_union.centroid
    print(f"\nModel centroid: {centroid.y:.4f}N, {centroid.x:.4f}W")

    # Download HUC12 watershed containing centroid
    wbd = WBD("huc12")
    huc12 = wbd.bygeom(Point(centroid.x, centroid.y), geo_crs="EPSG:4326")

    print(f"\nHUC12 Watershed:")
    print(f"  HUC ID: {huc12.iloc[0]['huc12']}")
    print(f"  Name: {huc12.iloc[0]['name']}")
    print(f"  Area: {huc12.iloc[0]['areasqkm']:.1f} sq km ({huc12.iloc[0]['areasqkm'] * 0.386102:.1f} sq mi)")

    HUC12_AVAILABLE = True

except ImportError:
    print("pygeohydro not available - skipping HUC12 analysis")
    print("Install with: pip install pygeohydro")
    HUC12_AVAILABLE = False
except Exception as e:
    print(f"Error downloading HUC12: {e}")
    HUC12_AVAILABLE = False
Text Only
pygeohydro not available - skipping HUC12 analysis
Install with: pip install pygeohydro
Python
# =============================================================================
# CALCULATE DRAINAGE COVERAGE
# =============================================================================
if HUC12_AVAILABLE and 'huc12' in dir():
    # Use equal-area projection for accurate area calculations
    # EPSG:5070 (Albers Equal Area Conic) is standard for US applications
    equal_area_crs = "EPSG:5070"

    # Project mesh areas to equal-area CRS for accurate area calculation
    mesh_areas_proj = mesh_areas.to_crs(equal_area_crs)
    mesh_total_area_sqkm = mesh_areas_proj.geometry.area.sum() / 1e6  # m^2 to km^2

    # Use reported area from WBD (already accurate from USGS)
    huc12_area_sqkm = huc12.iloc[0]['areasqkm']

    # Calculate coverage percentage
    coverage_pct = (mesh_total_area_sqkm / huc12_area_sqkm) * 100
    unmodeled_area_sqkm = huc12_area_sqkm - mesh_total_area_sqkm

    print("Drainage Coverage Analysis:")
    print(f"  2D Mesh Area: {mesh_total_area_sqkm:.2f} sq km ({mesh_total_area_sqkm * 0.386102:.2f} sq mi)")
    print(f"  HUC12 Area: {huc12_area_sqkm:.2f} sq km ({huc12_area_sqkm * 0.386102:.2f} sq mi)")
    print(f"  Coverage: {coverage_pct:.1f}%")
    print(f"  Unmodeled Area: {unmodeled_area_sqkm:.2f} sq km ({unmodeled_area_sqkm * 0.386102:.2f} sq mi)")

    if coverage_pct < 80:
        print(f"\nWARNING: Model covers only {coverage_pct:.1f}% of HUC12 drainage area")
        print("  Expect validation discrepancies due to unmodeled runoff contributions")
else:
    print("HUC12 not available - skipping coverage analysis")
Text Only
HUC12 not available - skipping coverage analysis
Python
# =============================================================================
# CREATE COVERAGE FIGURE
# =============================================================================
from matplotlib.lines import Line2D
from matplotlib.patches import Patch
from matplotlib.ticker import FuncFormatter

if HUC12_AVAILABLE and 'huc12' in dir():
    # Use projected CRS for proper visualization
    vis_crs = "EPSG:5070"  # Albers Equal Area Conic

    fig, ax = plt.subplots(figsize=(12, 10))

    # Project all data to visualization CRS
    huc12_proj = huc12.to_crs(vis_crs)
    mesh_areas_vis = mesh_areas.to_crs(vis_crs)
    legend_handles = [
        Patch(facecolor='lightblue', edgecolor='blue', alpha=0.3, label='HUC12 Watershed'),
        Patch(facecolor='green', edgecolor='darkgreen', alpha=0.5, label='2D Flow Areas'),
    ]
    gauge_roles = set()

    # Plot HUC12 boundary
    huc12_proj.plot(ax=ax, facecolor='lightblue', edgecolor='blue', linewidth=2, alpha=0.3)

    # Plot 2D mesh areas
    mesh_areas_vis.plot(ax=ax, facecolor='green', edgecolor='darkgreen', linewidth=1, alpha=0.5)

    # Add USGS gauge locations if we can retrieve them
    try:
        upstream_meta = get_gauge_metadata(UPSTREAM_GAUGE)
        validation_meta = get_gauge_metadata(VALIDATION_GAUGE)

        # Create gauge points in WGS84
        gauge_points = gpd.GeoDataFrame([
            {'site_id': UPSTREAM_GAUGE, 'name': upstream_meta['station_name'], 
             'geometry': Point(upstream_meta['longitude'], upstream_meta['latitude']), 'role': 'Upstream BC'},
            {'site_id': VALIDATION_GAUGE, 'name': validation_meta['station_name'],
             'geometry': Point(validation_meta['longitude'], validation_meta['latitude']), 'role': 'Validation'}
        ], crs="EPSG:4326")

        # Transform to visualization CRS
        gauge_points_proj = gauge_points.to_crs(vis_crs)

        # Plot gauges with different colors
        for idx, row in gauge_points_proj.iterrows():
            color = 'red' if row['role'] == 'Upstream BC' else 'orange'
            gauge_roles.add(row['role'])
            ax.scatter(row.geometry.x, row.geometry.y, c=color, s=150, marker='^', 
                      edgecolors='black', linewidths=1.5, zorder=5)
            ax.annotate(f"{row['site_id']}\n({row['role']})", 
                       (row.geometry.x, row.geometry.y), 
                       xytext=(10, 10), textcoords='offset points', fontsize=9,
                       bbox=dict(boxstyle='round,pad=0.2', facecolor='white', edgecolor='none', alpha=0.75))

    except Exception as e:
        print(f"Could not add gauge locations: {e}")

    if 'Upstream BC' in gauge_roles:
        legend_handles.append(Line2D([0], [0], marker='^', linestyle='None', markersize=10,
                                     markerfacecolor='red', markeredgecolor='black', label='Upstream BC Gauge'))
    if 'Validation' in gauge_roles:
        legend_handles.append(Line2D([0], [0], marker='^', linestyle='None', markersize=10,
                                     markerfacecolor='orange', markeredgecolor='black', label='Validation Gauge'))

    ax.set_title(f"Drainage Coverage Analysis\n{PROJECT_NAME} - HUC12: {huc12.iloc[0]['huc12']}", fontsize=14)
    ax.set_xlabel('Easting (km, EPSG:5070)')
    ax.set_ylabel('Northing (km, EPSG:5070)')
    ax.xaxis.set_major_formatter(FuncFormatter(lambda value, _: f'{value / 1000:,.0f}'))
    ax.yaxis.set_major_formatter(FuncFormatter(lambda value, _: f'{value / 1000:,.0f}'))
    ax.set_aspect('equal', adjustable='box')
    ax.legend(handles=legend_handles, loc='lower right', title='Map Layers', frameon=True, framealpha=0.95)

    # Add coverage annotation
    ax.annotate(f"Coverage: {coverage_pct:.1f}%\nUnmodeled: {unmodeled_area_sqkm:.1f} sq km",
               xy=(0.02, 0.98), xycoords='axes fraction', fontsize=10,
               verticalalignment='top', bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
    ax.annotate('N', xy=(0.08, 0.18), xytext=(0.08, 0.07), xycoords='axes fraction',
               ha='center', va='center', fontsize=11, fontweight='bold',
               arrowprops=dict(facecolor='black', edgecolor='black', width=2, headwidth=8))

    plt.tight_layout()
    plt.savefig(project_folder / "coverage_analysis.png", dpi=150)
    plt.show()

    print(f"\nCoverage figure saved to: {project_folder / 'coverage_analysis.png'}")
else:
    print("Skipping coverage figure (HUC12 not available)")
Text Only
Skipping coverage figure (HUC12 not available)
Python
# =============================================================================
# DOWNLOAD AORC PRECIPITATION DATA
# =============================================================================
# Create precipitation folder
precip_folder = project_folder / "Precipitation"
precip_folder.mkdir(exist_ok=True)

aorc_file = precip_folder / f"storm_{STORM_DATE}.nc"

print(f"Downloading AORC precipitation data...")
print(f"  Start: {SIM_START} 00:00")
print(f"  End: {SIM_END} 00:00")
print(f"  Output: {aorc_file}")

try:
    output_path = PrecipAorc.download(
        bounds=bounds,
        start_time=f"{SIM_START} 00:00",
        end_time=f"{SIM_END} 00:00",
        output_path=aorc_file,
        target_crs="EPSG:5070",  # SHG (Standard Hydrologic Grid) for HEC-RAS
        resolution=2000.0  # 2km resolution (standard SHG)
    )
    print(f"\nAORC data downloaded successfully: {output_path}")

    # Get file size
    file_size_mb = output_path.stat().st_size / (1024 * 1024)
    print(f"  File size: {file_size_mb:.2f} MB")

except Exception as e:
    print(f"\nError downloading AORC data: {e}")
    print("Continuing with workflow - manual precipitation setup may be needed")
    aorc_file = None
Text Only
Downloading AORC precipitation data...
  Start: 2020-12-22 00:00
  End: 2020-12-27 00:00
  Output: C:\Users\bill\.config\superpowers\worktrees\ras-commander-pr-docs\codex-pr251-timeseries-alignment\examples\example_projects\BaldEagleCrkMulti2D_914_historical\Precipitation\storm_20201224.nc


2026-07-08 21:38:38 - ras_commander.precip.PrecipAorc - INFO - Downloading AORC APCP_surface: 2020-12-22 00:00:00 to 2020-12-27 00:00:00, bounds W=-77.7085, S=41.0102, E=-77.2511, N=41.2197


2026-07-08 21:38:44 - ras_commander.precip.PrecipAorc - INFO - AORC download complete: storm_20201224.nc (0.1 MB, 144 timesteps)



AORC data downloaded successfully: C:\Users\bill\.config\superpowers\worktrees\ras-commander-pr-docs\codex-pr251-timeseries-alignment\examples\example_projects\BaldEagleCrkMulti2D_914_historical\Precipitation\storm_20201224.nc
  File size: 0.12 MB
Python
# =============================================================================
# RETRIEVE USGS UPSTREAM BOUNDARY CONDITION DATA
# =============================================================================
print(f"Retrieving upstream BC flow data from gauge {UPSTREAM_GAUGE}...")

# Get gauge metadata
upstream_meta = get_gauge_metadata(UPSTREAM_GAUGE)
print(f"  Station: {upstream_meta['station_name']}")
print(f"  Drainage area: {upstream_meta.get('drainage_area_sqmi', 'N/A')} sq mi")

# Retrieve flow data for simulation period
upstream_flow = retrieve_flow_data(
    site_id=UPSTREAM_GAUGE,
    start_datetime=SIM_START,
    end_datetime=SIM_END,
    data_type='iv'  # Instantaneous values
)

# Remove timezone to match HEC-RAS (timezone-naive)
upstream_flow['datetime'] = pd.to_datetime(upstream_flow['datetime']).dt.tz_localize(None)

print(f"\nRetrieved {len(upstream_flow)} records")
print(f"  Period: {upstream_flow['datetime'].min()} to {upstream_flow['datetime'].max()}")
print(f"  Flow range: {upstream_flow['value'].min():.0f} to {upstream_flow['value'].max():.0f} cfs")
Text Only
Retrieving upstream BC flow data from gauge 01547500...


2026-07-08 21:38:45 - httpx - INFO - HTTP Request: GET https://waterservices.usgs.gov/nwis/site?sites=01547500&siteOutput=Expanded&format=rdb "HTTP/1.1 200 "


  Station: Bald Eagle Creek at Blanchard, PA
  Drainage area: 339.0 sq mi


2026-07-08 21:38:45 - httpx - INFO - HTTP Request: GET https://waterservices.usgs.gov/nwis/iv?format=json&parameterCd=00060&startDT=2020-12-22&endDT=2020-12-27&sites=01547500 "HTTP/1.1 301 "


2026-07-08 21:38:46 - httpx - INFO - HTTP Request: GET https://nwis.waterservices.usgs.gov/nwis/iv/?format=json&parameterCd=00060&startDT=2020-12-22&endDT=2020-12-27&sites=01547500 "HTTP/1.1 200 "



Retrieved 576 records
  Period: 2020-12-22 05:00:00 to 2020-12-28 04:45:00
  Flow range: 1 to 2040 cfs
Python
# =============================================================================
# RETRIEVE USGS VALIDATION DATA
# =============================================================================
print(f"Retrieving validation data from gauge {VALIDATION_GAUGE}...")

# Get gauge metadata
try:
    validation_meta = get_gauge_metadata(VALIDATION_GAUGE)
    print(f"  Station: {validation_meta['station_name']}")
    print(f"  Drainage area: {validation_meta.get('drainage_area_sqmi', 'N/A')} sq mi")
except Exception as e:
    print(f"  Could not retrieve metadata: {e}")
    validation_meta = {'station_name': f'USGS {VALIDATION_GAUGE}'}

# Retrieve flow data for validation
try:
    validation_flow = retrieve_flow_data(
        site_id=VALIDATION_GAUGE,
        start_datetime=SIM_START,
        end_datetime=SIM_END,
        data_type='iv'
    )

    # Remove timezone
    validation_flow['datetime'] = pd.to_datetime(validation_flow['datetime']).dt.tz_localize(None)

    print(f"\nRetrieved {len(validation_flow)} flow records")
    print(f"  Period: {validation_flow['datetime'].min()} to {validation_flow['datetime'].max()}")
    print(f"  Flow range: {validation_flow['value'].min():.0f} to {validation_flow['value'].max():.0f} cfs")

except Exception as e:
    print(f"Error retrieving validation flow data: {e}")
    validation_flow = None

# Retrieve stage data for validation
try:
    validation_stage = retrieve_stage_data(
        site_id=VALIDATION_GAUGE,
        start_datetime=SIM_START,
        end_datetime=SIM_END,
        data_type='iv'
    )

    # Remove timezone
    validation_stage['datetime'] = pd.to_datetime(validation_stage['datetime']).dt.tz_localize(None)

    print(f"\nRetrieved {len(validation_stage)} stage records")
    print(f"  Stage range: {validation_stage['value'].min():.2f} to {validation_stage['value'].max():.2f} ft")

except Exception as e:
    print(f"Stage data not available for gauge {VALIDATION_GAUGE}: {e}")
    validation_stage = None
Text Only
Retrieving validation data from gauge 01548005...


2026-07-08 21:38:46 - httpx - INFO - HTTP Request: GET https://waterservices.usgs.gov/nwis/site?sites=01548005&siteOutput=Expanded&format=rdb "HTTP/1.1 200 "


2026-07-08 21:38:47 - httpx - INFO - HTTP Request: GET https://waterservices.usgs.gov/nwis/iv?format=json&parameterCd=00060&startDT=2020-12-22&endDT=2020-12-27&sites=01548005 "HTTP/1.1 301 "


  Station: Bald Eagle Creek near Beech Creek Station, PA
  Drainage area: 562.0 sq mi


2026-07-08 21:38:47 - httpx - INFO - HTTP Request: GET https://nwis.waterservices.usgs.gov/nwis/iv/?format=json&parameterCd=00060&startDT=2020-12-22&endDT=2020-12-27&sites=01548005 "HTTP/1.1 200 "



Retrieved 576 flow records
  Period: 2020-12-22 05:00:00 to 2020-12-28 04:45:00
  Flow range: 378 to 4340 cfs


2026-07-08 21:38:47 - httpx - INFO - HTTP Request: GET https://waterservices.usgs.gov/nwis/iv?format=json&parameterCd=00065&startDT=2020-12-22&endDT=2020-12-27&sites=01548005 "HTTP/1.1 301 "


2026-07-08 21:38:48 - httpx - INFO - HTTP Request: GET https://nwis.waterservices.usgs.gov/nwis/iv/?format=json&parameterCd=00065&startDT=2020-12-22&endDT=2020-12-27&sites=01548005 "HTTP/1.1 200 "



Retrieved 576 stage records
  Stage range: 6.67 to 11.78 ft
Python
# =============================================================================
# PLOT UPSTREAM BC HYDROGRAPH
# =============================================================================
fig, ax = plt.subplots(figsize=(12, 5))

ax.plot(upstream_flow['datetime'], upstream_flow['value'], 'b-', linewidth=1, label='USGS Flow')
ax.fill_between(upstream_flow['datetime'], 0, upstream_flow['value'], alpha=0.2)

# Highlight storm period
storm_start = pd.Timestamp('2020-12-24')
storm_end = pd.Timestamp('2020-12-26')
ax.axvspan(storm_start, storm_end, alpha=0.1, color='red', label='Storm Period')

ax.set_xlabel('Date')
ax.set_ylabel('Flow (cfs)')
ax.set_title(f'Upstream Boundary Condition: USGS {UPSTREAM_GAUGE}\n{upstream_meta["station_name"]}')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Calculate peak flow during storm
storm_data = upstream_flow[(upstream_flow['datetime'] >= storm_start) & (upstream_flow['datetime'] <= storm_end)]
print(f"\nStorm Period Statistics ({storm_start.strftime('%Y-%m-%d')} to {storm_end.strftime('%Y-%m-%d')})")
print(f"  Peak flow: {storm_data['value'].max():.0f} cfs")
print(f"  Peak time: {storm_data.loc[storm_data['value'].idxmax(), 'datetime']}")

png

Text Only
Storm Period Statistics (2020-12-24 to 2020-12-26)
  Peak flow: 544 cfs
  Peak time: 2020-12-24 21:30:00
Python
# =============================================================================
# CLONE PLAN AND UNSTEADY FILE FOR THE 2020 EVENT
# =============================================================================
print(f"Cloning template plan {TEMPLATE_PLAN} and its unsteady file for {STORM_DATE}...")

sim_start_dt = pd.to_datetime(SIM_START).to_pydatetime()
sim_end_dt = pd.to_datetime(SIM_END).to_pydatetime()
precip_rel_path = Path("Precipitation") / f"storm_{STORM_DATE}.nc"
precip_full_path = project_folder / precip_rel_path

if not precip_full_path.exists():
    raise FileNotFoundError(f"AORC NetCDF file not found: {precip_full_path}")

template_plan_row = ras.plan_df[
    ras.plan_df["plan_number"].astype(str).str.zfill(2).eq(str(TEMPLATE_PLAN).zfill(2))
].iloc[0]
template_unsteady = str(template_plan_row["unsteady_number"]).zfill(2)

# Keep the template 2018 unsteady file intact and isolate the 2020 event in a clone.
new_unsteady = RasPlan.clone_unsteady(
    template_unsteady,
    new_title=f"AORC {STORM_DATE}",
    ras_object=ras,
)
new_plan = RasPlan.clone_plan(
    TEMPLATE_PLAN,
    new_plan_shortid="storm12",
    new_title="AORC Dec 2020 Validation",
    unsteady_flow=new_unsteady,
    ras_object=ras,
)

RasPlan.update_simulation_date(
    new_plan,
    start_date=sim_start_dt,
    end_date=sim_end_dt,
    ras_object=ras,
)
RasUnsteady.set_gridded_precipitation(
    unsteady_file=new_unsteady,
    netcdf_path=precip_rel_path,
    interpolation="Nearest",
    ras_object=ras,
)

# Re-initialize to refresh plan, unsteady, and boundary DataFrames after edits.
ras = init_ras_project(project_folder, RAS_VERSION)

print(f"Created plan {new_plan} using cloned unsteady file u{new_unsteady}")
print(f"Simulation window: {SIM_START} 00:00 to {SIM_END} 00:00")
print(f"AORC precipitation source: {precip_rel_path}")

plan_cols = ["plan_number", "Plan Title", "Simulation Date", "unsteady_number"]
print("\nUpdated event plan:")
display(ras.plan_df.loc[ras.plan_df["plan_number"].astype(str).str.zfill(2).eq(str(new_plan).zfill(2)), plan_cols])
Text Only
Cloning template plan 06 and its unsteady file for 20201224...


2026-07-08 21:38:49 - ras_commander.RasPlan - INFO - Updated simulation date in plan file: BaldEagleDamBrk.p07


2026-07-08 21:38:50 - ras_commander.RasUnsteady - INFO - Configured gridded precipitation in BaldEagleDamBrk.u04: source=.\Precipitation\storm_20201224.nc, interpolation=Nearest, dataset=APCP_surface


2026-07-08 21:38:50 - ras_commander.RasUnsteady - INFO - Imported gridded precipitation into BaldEagleDamBrk.u04.hdf: 144 timesteps, 352 cells, range=0.0-80.3 mm


2026-07-08 21:38:50 - ras_commander.RasPrj - INFO - ras-commander v0.98.2 | An open-source project of CLB Engineering Corporation (https://clbengineering.com/) | Docs: https://rascommander.info | GitHub: https://github.com/gpt-cmdr/ras-commander


2026-07-08 21:38:50 - ras_commander.RasPrj - INFO - Project initialized: BaldEagleDamBrk


2026-07-08 21:38:50 - ras_commander.RasPrj - INFO - Using HEC-RAS version 7.0


2026-07-08 21:38:50 - ras_commander.RasPrj - INFO - 
═══════════════════════════════════════════════════════════════════════
ras-commander | HEC-RAS Automation Library
Docs: https://rascommander.info/
Repo: https://github.com/gpt-cmdr/ras-commander
LLM agents: https://rascommander.info/llms.txt
═══════════════════════════════════════════════════════════════════════

PROJECT DATAFRAMES (single source of truth — use these, not file globbing):
  ras.plan_df        Plans, HDF paths, geometry/flow associations
  ras.geom_df        Geometry files and HDF preprocessor paths
  ras.flow_df        Steady flow files
  ras.unsteady_df    Unsteady flow files and configurations
  ras.boundaries_df  Boundary conditions (type, name, location)
  ras.results_df     Lightweight HDF results summaries
  ras.rasmap_df      RASMapper layers, terrain, land cover paths

KEY APIS (static classes — call directly, never instantiate):
  Execution:    RasCmdr.compute_plan() / compute_parallel() / compute_test_mode()
  Plan Files:   RasPlan.clone_plan() / clone_geom() / set_geom()
  Unsteady:     RasUnsteady — IC/BC management, gate openings, precipitation
  Geometry:     GeomCrossSection, GeomBridge, GeomStorage, GeomLateral, GeomMesh
  HDF Results:  HdfResultsPlan.get_wse() / get_compute_messages()
                HdfResultsMesh.get_mesh_max_ws() / get_mesh_cells_timeseries()
                HdfMesh.get_mesh_cell_points()
  QA/QC:        RasCheck.run_check() / RasFixit (geometry repair)
  DSS:          RasDss.get_timeseries() / check_pathname()
  USGS:         UsgsGaugeSpatial, GaugeMatcher, RasUsgsBoundaryGeneration
  Precipitation: StormGenerator, Atlas14Storm, PrecipAorc, Atlas14Variance
  Terrain:      RasTerrain.create_terrain_hdf() / RasTerrainMod

MULTI-PROJECT: Pass ras_object= to all API calls when using local RasPrj instances.

EXAMPLES: 100+ notebooks in examples/ (100s=execution, 200s=geometry, 300s=unsteady,
  400s=HDF results, 500s=remote, 800s=QA/QC, 900s=data integration).
  Review relevant notebooks before assembling new workflows.

PLATFORM: Most HEC-RAS operations require Windows. Linux/Wine support for
  headless execution, data access, geometry modification, and preprocessing
  is available via RasProcess (HEC-RAS 6.6+). See ras_commander/RasProcess.py.
  Remote distributed execution: ras_commander/remote/ (PsExec, Docker, SSH, cloud).
═══════════════════════════════════════════════════════════════════════


Created plan 07 using cloned unsteady file u04
Simulation window: 2020-12-22 00:00 to 2020-12-27 00:00
AORC precipitation source: Precipitation\storm_20201224.nc

Updated event plan:
plan_number Plan Title Simulation Date unsteady_number
11 07 AORC Dec 2020 Validation 22DEC2020,0000,27DEC2020,0000 04
Python
# =============================================================================
# VERIFY EVENT INPUT ALIGNMENT
# =============================================================================
plan_row = ras.plan_df.loc[
    ras.plan_df["plan_number"].astype(str).str.zfill(2).eq(str(new_plan).zfill(2))
].iloc[0]
active_unsteady = str(plan_row["unsteady_number"]).zfill(2)

print("Event input alignment check")
print("=" * 70)
print(f"Plan {new_plan} simulation date: {plan_row['Simulation Date']}")
print(f"Plan {new_plan} unsteady file: u{active_unsteady}")
print(f"Precipitation NetCDF: {precip_rel_path}")
print(f"Precipitation file exists: {precip_full_path.exists()}")

boundary_cols = [
    col for col in [
        "unsteady_number", "bc_type", "river_reach_name", "river_station",
        "SA/2D Connection", "Gate Name", "Interval", "Use DSS", "DSS Path"
    ]
    if col in ras.boundaries_df.columns
]
active_boundaries = ras.boundaries_df[
    ras.boundaries_df["unsteady_number"].astype(str).str.zfill(2).eq(active_unsteady)
].copy()
print(f"\nActive boundary rows for u{active_unsteady}: {len(active_boundaries)}")
display(active_boundaries[boundary_cols].head(12) if boundary_cols else active_boundaries.head(12))

print("\nConstraints and assumptions:")
print("- AORC NetCDF precipitation is the active event forcing for this 2020 run.")
print("- The cloned plan, cloned unsteady file, and precipitation source are aligned to the same event window.")
print("- USGS records are used for validation context in this notebook; replacing model boundary hydrographs with observed USGS data is a separate calibration/validation step.")
print("- Gate/control operations are preserved from the template model unless explicitly replaced with observed operations data.")
Text Only
Event input alignment check
======================================================================
Plan 07 simulation date: 22DEC2020,0000,27DEC2020,0000
Plan 07 unsteady file: u04
Precipitation NetCDF: Precipitation\storm_20201224.nc
Precipitation file exists: True

Active boundary rows for u04: 4
unsteady_number bc_type river_reach_name river_station Interval Use DSS DSS Path
51 04 Normal Depth NaN NaN NaN
52 04 Flow Hydrograph 1HOUR False
53 04 Normal Depth NaN NaN NaN
54 04 Gate Opening NaN NaN NaN
Text Only
Constraints and assumptions:
- AORC NetCDF precipitation is the active event forcing for this 2020 run.
- The cloned plan, cloned unsteady file, and precipitation source are aligned to the same event window.
- USGS records are used for validation context in this notebook; replacing model boundary hydrographs with observed USGS data is a separate calibration/validation step.
- Gate/control operations are preserved from the template model unless explicitly replaced with observed operations data.

Event Forcing Constraints and Assumptions

This notebook now configures the cloned HEC-RAS plan programmatically before execution:

  • The 2020 validation plan is isolated in a cloned unsteady-flow file so the original 2018 template event remains unchanged.
  • The plan simulation window is 2020-12-22 00:00 through 2020-12-27 00:00.
  • The active precipitation source is the downloaded AORC NetCDF file, Precipitation/storm_20201224.nc.
  • Existing model boundary hydrographs and gate/control operations are preserved from the template model unless explicitly replaced with observed data.
  • USGS records retrieved above are used for validation context. Injecting observed upstream/downstream hydrographs into the active unsteady file is a separate boundary-condition calibration step and is not assumed here.
Python
# =============================================================================
# RUN MODEL (Skip if manual configuration needed)
# =============================================================================
RUN_MODEL = True  # Set to True after manual configuration

if RUN_MODEL:
    print(f"Running plan {new_plan}...")

    try:
        RasCmdr.compute_plan(
            plan_number=new_plan,
            num_cores=4,
            force_rerun=True,
            verify=True,
            ras_object=ras
        )
        print("Model execution complete")

        # Re-initialize to pick up results
        ras = init_ras_project(project_folder, RAS_VERSION)

    except Exception as e:
        print(f"Model execution failed: {e}")
        print("\nCheck that:")
        print("  1. The cloned unsteady file points to the 2020 AORC NetCDF")
        print("  2. Simulation dates match the precipitation data window")
        print("  3. Boundary and gate/control time series cover the plan period")
else:
    print("Model execution skipped - set RUN_MODEL = True after configuration")
    print("\nAlternatively, use existing results from pre-run plan...")
Text Only
Running plan 07...


Model execution complete


2026-07-08 21:44:07 - ras_commander.RasPrj - INFO - ras-commander v0.98.2 | An open-source project of CLB Engineering Corporation (https://clbengineering.com/) | Docs: https://rascommander.info | GitHub: https://github.com/gpt-cmdr/ras-commander


2026-07-08 21:44:07 - ras_commander.RasPrj - INFO - Project initialized: BaldEagleDamBrk


2026-07-08 21:44:07 - ras_commander.RasPrj - INFO - Using HEC-RAS version 7.0


2026-07-08 21:44:07 - ras_commander.RasPrj - INFO - 
═══════════════════════════════════════════════════════════════════════
ras-commander | HEC-RAS Automation Library
Docs: https://rascommander.info/
Repo: https://github.com/gpt-cmdr/ras-commander
LLM agents: https://rascommander.info/llms.txt
═══════════════════════════════════════════════════════════════════════

PROJECT DATAFRAMES (single source of truth — use these, not file globbing):
  ras.plan_df        Plans, HDF paths, geometry/flow associations
  ras.geom_df        Geometry files and HDF preprocessor paths
  ras.flow_df        Steady flow files
  ras.unsteady_df    Unsteady flow files and configurations
  ras.boundaries_df  Boundary conditions (type, name, location)
  ras.results_df     Lightweight HDF results summaries
  ras.rasmap_df      RASMapper layers, terrain, land cover paths

KEY APIS (static classes — call directly, never instantiate):
  Execution:    RasCmdr.compute_plan() / compute_parallel() / compute_test_mode()
  Plan Files:   RasPlan.clone_plan() / clone_geom() / set_geom()
  Unsteady:     RasUnsteady — IC/BC management, gate openings, precipitation
  Geometry:     GeomCrossSection, GeomBridge, GeomStorage, GeomLateral, GeomMesh
  HDF Results:  HdfResultsPlan.get_wse() / get_compute_messages()
                HdfResultsMesh.get_mesh_max_ws() / get_mesh_cells_timeseries()
                HdfMesh.get_mesh_cell_points()
  QA/QC:        RasCheck.run_check() / RasFixit (geometry repair)
  DSS:          RasDss.get_timeseries() / check_pathname()
  USGS:         UsgsGaugeSpatial, GaugeMatcher, RasUsgsBoundaryGeneration
  Precipitation: StormGenerator, Atlas14Storm, PrecipAorc, Atlas14Variance
  Terrain:      RasTerrain.create_terrain_hdf() / RasTerrainMod

MULTI-PROJECT: Pass ras_object= to all API calls when using local RasPrj instances.

EXAMPLES: 100+ notebooks in examples/ (100s=execution, 200s=geometry, 300s=unsteady,
  400s=HDF results, 500s=remote, 800s=QA/QC, 900s=data integration).
  Review relevant notebooks before assembling new workflows.

PLATFORM: Most HEC-RAS operations require Windows. Linux/Wine support for
  headless execution, data access, geometry modification, and preprocessing
  is available via RasProcess (HEC-RAS 6.6+). See ras_commander/RasProcess.py.
  Remote distributed execution: ras_commander/remote/ (PsExec, Docker, SSH, cloud).
═══════════════════════════════════════════════════════════════════════
Python
# =============================================================================
# USE EXISTING PLAN RESULTS (If available)
# =============================================================================
# Check if any plan has HDF results we can use for demonstration
hdf_path = None

for idx, row in ras.plan_df.iterrows():
    hdf_path_str = row.get('HDF_Results_Path')
    if isinstance(hdf_path_str, str) and Path(hdf_path_str).exists():
        hdf_path = Path(hdf_path_str)
        plan_number = row['plan_number']
        print(f"Found existing results: Plan {plan_number}")
        print(f"  HDF: {hdf_path.name}")
        break

if hdf_path is None:
    print("No existing HDF results found.")
    print("Run the model first or use an existing plan with results.")
Text Only
Found existing results: Plan 07
  HDF: BaldEagleDamBrk.p07.hdf

=============================================================================

VALIDATION: EXTRACT AND COMPARE MODELED VS OBSERVED DATA

=============================================================================

The cells below extract modeled results from the HDF file and compare them

against observed USGS data. This section requires:

1. A completed HEC-RAS run with HDF results (check cell above)

2. validation_meta defined (from cell 11)

3. validation_flow/validation_stage defined (from cell 11)

Python
# =============================================================================
# EXTRACT MODELED RESULTS
# =============================================================================
import h5py

modeled_flow_df = None
modeled_stage_df = None
modeled_df = None  # retained for compatibility with older narrative text
stage_extraction_meta = {}

XSEC_OUTPUT_PATH = (
    "Results/Unsteady/Output/Output Blocks/Base Output/"
    "Unsteady Time Series/Cross Sections"
)
REF_LINE_OUTPUT_PATH = (
    "Results/Unsteady/Output/Output Blocks/Base Output/"
    "Unsteady Time Series/Reference Lines"
)
MESH_TS_ROOT = (
    "Results/Unsteady/Output/Output Blocks/Base Output/"
    "Unsteady Time Series/2D Flow Areas"
)
TIME_STAMP_PATH = (
    "Results/Unsteady/Output/Output Blocks/Base Output/"
    "Unsteady Time Series/Time Date Stamp"
)
PROJECT_CRS = "EPSG:2271"  # Bald Eagle Creek example project CRS (PA StatePlane North, US feet)


def _decode_ras_timestamps(raw_times):
    """Decode HEC-RAS timestamp byte strings into pandas datetimes."""
    decoded = [
        t.decode("utf-8").strip() if isinstance(t, (bytes, bytearray)) else str(t).strip()
        for t in raw_times
    ]
    return pd.to_datetime(decoded, format="%d%b%Y %H:%M:%S")


def _extract_nearest_cell_wse(plan_hdf_path, gauge_meta, project_crs=PROJECT_CRS):
    """Extract WSE time series from the 2D mesh cell nearest the validation gauge."""
    required_meta = {"latitude", "longitude"}
    missing_meta = required_meta.difference(gauge_meta or {})
    if missing_meta:
        raise ValueError(f"Validation gauge metadata is missing: {sorted(missing_meta)}")

    with h5py.File(plan_hdf_path, "r") as hdf:
        if MESH_TS_ROOT not in hdf:
            raise KeyError(f"2D mesh time-series group not found: {MESH_TS_ROOT}")
        if TIME_STAMP_PATH not in hdf:
            raise KeyError(f"Time stamp dataset not found: {TIME_STAMP_PATH}")

        mesh_names = [
            mesh_name for mesh_name in hdf[MESH_TS_ROOT].keys()
            if f"{MESH_TS_ROOT}/{mesh_name}/Water Surface" in hdf
        ]
        if not mesh_names:
            raise KeyError("No 2D mesh Water Surface time-series datasets found")

        mesh_name = mesh_names[0]
        centers_path = f"Geometry/2D Flow Areas/{mesh_name}/Cells Center Coordinate"
        wse_path = f"{MESH_TS_ROOT}/{mesh_name}/Water Surface"
        if centers_path not in hdf:
            raise KeyError(f"Cell-center coordinate dataset not found: {centers_path}")

        centers = hdf[centers_path][:]
        gauge_wgs84 = gpd.GeoDataFrame(
            [{"site_id": VALIDATION_GAUGE}],
            geometry=[Point(float(gauge_meta["longitude"]), float(gauge_meta["latitude"]))],
            crs="EPSG:4326",
        )
        gauge_projected = gauge_wgs84.to_crs(project_crs).geometry.iloc[0]
        distances = np.hypot(centers[:, 0] - gauge_projected.x, centers[:, 1] - gauge_projected.y)
        nearest_cell = int(np.nanargmin(distances))
        distance_ft = float(distances[nearest_cell])
        distance_m = distance_ft * 0.3048006096

        wse_values = hdf[wse_path][:, nearest_cell].astype(float)
        wse_values = np.where(wse_values < -1.0e20, np.nan, wse_values)
        times = _decode_ras_timestamps(hdf[TIME_STAMP_PATH][:])

    modeled_stage = pd.DataFrame({
        "datetime": times,
        "value": wse_values,
    }).dropna()

    meta = {
        "mesh_name": mesh_name,
        "cell_id": nearest_cell,
        "distance_m": distance_m,
        "distance_ft": distance_ft,
        "project_crs": project_crs,
    }
    return modeled_stage, meta


if hdf_path and hdf_path.exists():
    print(f"Extracting modeled results from {hdf_path.name}...")

    # Try 1D cross-section flow only if that output group exists.
    with h5py.File(hdf_path, "r") as hdf:
        has_xsec_output = XSEC_OUTPUT_PATH in hdf
        has_ref_line_output = REF_LINE_OUTPUT_PATH in hdf

    if has_xsec_output:
        print("Attempting to extract 1D cross-section flow results...")
        xs_data = HdfResultsXsec.get_xsec_timeseries(hdf_path)
        xs_names = xs_data.coords["cross_section"].values.tolist()
        station_values = xs_data.coords["Station"].values.tolist()
        time_values = xs_data.coords["time"].values

        target_station = float(station_values[-1]) + (float(station_values[0]) - float(station_values[-1])) * 0.3
        closest_idx = min(range(len(station_values)), key=lambda i: abs(float(station_values[i]) - target_station))
        validation_xs = xs_names[closest_idx]
        modeled_flow_ts = xs_data["Flow"].sel(cross_section=validation_xs)

        modeled_flow_df = pd.DataFrame({
            "datetime": pd.to_datetime(time_values),
            "value": modeled_flow_ts.values,
        })
        modeled_df = modeled_flow_df

        print("\nModeled flow results (1D cross-section):")
        print(f"  Cross sections: {len(xs_names)}")
        print(f"  Validation cross section: {validation_xs} (station {station_values[closest_idx]})")
        print(f"  Period: {modeled_flow_df['datetime'].min()} to {modeled_flow_df['datetime'].max()}")
        print(f"  Flow range: {modeled_flow_df['value'].min():.1f} to {modeled_flow_df['value'].max():.1f} cfs")
    else:
        print("1D cross-section flow output is not present in this plan HDF.")

    # Try 2D reference-line flow next. Empty reference-line output is expected in this example project.
    if modeled_flow_df is None and has_ref_line_output:
        print("\nAttempting to extract 2D reference-line flow results...")
        ref_lines_data = HdfResultsXsec.get_ref_lines_timeseries(hdf_path)
        if ref_lines_data and len(ref_lines_data.data_vars) > 0 and "Flow" in ref_lines_data.data_vars:
            ref_line_names = ref_lines_data.coords["ref_line"].values.tolist()
            validation_ref_line = ref_line_names[0]
            modeled_flow_ts = ref_lines_data["Flow"].sel(ref_line=validation_ref_line)
            modeled_flow_df = pd.DataFrame({
                "datetime": pd.to_datetime(ref_lines_data.coords["time"].values),
                "value": modeled_flow_ts.values,
            })
            modeled_df = modeled_flow_df

            print("\nModeled flow results (2D reference line):")
            print(f"  Reference line: {validation_ref_line}")
            print(f"  Period: {modeled_flow_df['datetime'].min()} to {modeled_flow_df['datetime'].max()}")
            print(f"  Flow range: {modeled_flow_df['value'].min():.1f} to {modeled_flow_df['value'].max():.1f} cfs")
        else:
            print("2D reference-line flow output is empty or does not include Flow.")
    elif modeled_flow_df is None:
        print("2D reference-line flow output is not present in this plan HDF.")

    # If no defensible flow time series exists, extract nearest-cell WSE for stage validation.
    if modeled_flow_df is None:
        print("\nFlow validation requires a cross section, 2D reference line, or face-flow integration.")
        print("Falling back to nearest-cell 2D Water Surface Elevation for stage validation...")
        try:
            modeled_stage_df, stage_extraction_meta = _extract_nearest_cell_wse(hdf_path, validation_meta)
            modeled_df = modeled_stage_df
            print("\nModeled stage results (nearest 2D mesh cell WSE):")
            print(f"  Mesh area: {stage_extraction_meta['mesh_name']}")
            print(f"  Cell ID: {stage_extraction_meta['cell_id']}")
            print(
                f"  Distance from validation gauge: "
                f"{stage_extraction_meta['distance_ft']:.0f} ft "
                f"({stage_extraction_meta['distance_m']:.1f} m)"
            )
            print(f"  Period: {modeled_stage_df['datetime'].min()} to {modeled_stage_df['datetime'].max()}")
            print(f"  WSE range: {modeled_stage_df['value'].min():.2f} to {modeled_stage_df['value'].max():.2f} ft")
        except Exception as stage_error:
            print(f"  Could not extract nearest-cell stage series: {stage_error}")
            modeled_stage_df = None

    if modeled_flow_df is None and modeled_stage_df is None:
        print("\n" + "="*70)
        print("WARNING: Could not extract modeled flow or stage time series from the model results.")
        print("="*70)
        print("Options:")
        print("  1. Add a reference line in HEC-RAS at the validation gauge location")
        print("  2. Add a reference point or extract a reviewed 2D mesh-cell WSE location")
        print("  3. Calculate flow from face velocities with a reviewed cutline")
        print("="*70)
else:
    print("No HDF results available for extraction")
Text Only
Extracting modeled results from BaldEagleDamBrk.p07.hdf...
1D cross-section flow output is not present in this plan HDF.
2D reference-line flow output is not present in this plan HDF.

Flow validation requires a cross section, 2D reference line, or face-flow integration.
Falling back to nearest-cell 2D Water Surface Elevation for stage validation...



Modeled stage results (nearest 2D mesh cell WSE):
  Mesh area: BaldEagleCr
  Cell ID: 10562
  Distance from validation gauge: 97 ft (29.6 m)
  Period: 2020-12-22 00:00:00 to 2020-12-27 00:00:00
  WSE range: 562.58 to 573.83 ft
Python
# =============================================================================
# ALIGN AND CALCULATE VALIDATION METRICS
# =============================================================================
aligned_flow = None
flow_metrics = None
aligned_stage = None
stage_metrics = None
observed_stage_wse = None


def _align_to_model_times(modeled_series, observed_series):
    """Interpolate observed values to modeled timestamps over their overlap."""
    model = modeled_series[["datetime", "value"]].dropna().sort_values("datetime").copy()
    obs = observed_series[["datetime", "value"]].dropna().sort_values("datetime").copy()
    model["datetime"] = pd.to_datetime(model["datetime"])
    obs["datetime"] = pd.to_datetime(obs["datetime"])

    overlap_start = max(model["datetime"].min(), obs["datetime"].min())
    overlap_end = min(model["datetime"].max(), obs["datetime"].max())
    if overlap_start >= overlap_end:
        raise ValueError(
            f"No overlapping time period. Model: {model['datetime'].min()} to {model['datetime'].max()}, "
            f"Observed: {obs['datetime'].min()} to {obs['datetime'].max()}"
        )

    model = model[(model["datetime"] >= overlap_start) & (model["datetime"] <= overlap_end)].copy()
    obs = obs[(obs["datetime"] >= overlap_start) & (obs["datetime"] <= overlap_end)].copy()

    model_seconds = model["datetime"].astype("int64") / 1e9
    obs_seconds = obs["datetime"].astype("int64") / 1e9
    observed_interp = np.interp(model_seconds, obs_seconds, obs["value"].astype(float))

    return pd.DataFrame({
        "datetime": model["datetime"].values,
        "modeled": model["value"].astype(float).values,
        "observed": observed_interp,
    })


if modeled_flow_df is not None and validation_flow is not None:
    print("Aligning modeled and observed flow timeseries...")

    aligned_flow = align_timeseries(
        modeled_df=modeled_flow_df,
        observed_df=validation_flow,
    )

    print(f"\nAligned {len(aligned_flow)} timesteps")
    print(f"  Period: {aligned_flow['datetime'].min()} to {aligned_flow['datetime'].max()}")
    print(f"  Modeled range: {aligned_flow['modeled'].min():.1f} to {aligned_flow['modeled'].max():.1f} cfs")
    print(f"  Observed range: {aligned_flow['observed'].min():.1f} to {aligned_flow['observed'].max():.1f} cfs")

    print("\nCalculating flow validation metrics...")
    flow_metrics = calculate_all_metrics(
        observed=aligned_flow["observed"],
        modeled=aligned_flow["modeled"],
        time_index=aligned_flow["datetime"],
    )

    print("\n" + "="*60)
    print("FLOW VALIDATION METRICS")
    print("="*60)
    print(f"Nash-Sutcliffe Efficiency (NSE): {flow_metrics['nse']:.3f}")
    print(f"Kling-Gupta Efficiency (KGE): {flow_metrics['kge']:.3f}")
    print(f"Peak Flow Error: {flow_metrics['peak_error_pct']:.1f}%")
    print(f"Volume Error: {flow_metrics['vol_error_pct']:.1f}%")
    print(f"RMSE: {flow_metrics['rmse']:.1f} cfs")
elif modeled_stage_df is not None and validation_stage is not None:
    print("Aligning modeled WSE and observed USGS stage timeseries...")

    gage_datum_ft = validation_meta.get("gage_datum_ft") if validation_meta else None
    if gage_datum_ft is None:
        print("Cannot calculate stage metrics because the USGS gage datum is unavailable.")
    else:
        observed_stage_wse = validation_stage.copy()
        observed_stage_wse["value"] = pd.to_numeric(observed_stage_wse["value"], errors="coerce") + float(gage_datum_ft)

        aligned_stage = _align_to_model_times(modeled_stage_df, observed_stage_wse)
        residuals = aligned_stage["modeled"] - aligned_stage["observed"]
        model_peak_idx = aligned_stage["modeled"].idxmax()
        obs_peak_idx = aligned_stage["observed"].idxmax()
        peak_time_error_hr = (
            aligned_stage.loc[model_peak_idx, "datetime"] - aligned_stage.loc[obs_peak_idx, "datetime"]
        ) / pd.Timedelta(hours=1)

        stage_metrics = {
            "gage_datum_ft": float(gage_datum_ft),
            "bias_ft": float(residuals.mean()),
            "rmse_ft": float(np.sqrt(np.mean(residuals ** 2))),
            "mae_ft": float(np.mean(np.abs(residuals))),
            "peak_error_ft": float(aligned_stage["modeled"].max() - aligned_stage["observed"].max()),
            "peak_time_error_hr": float(peak_time_error_hr),
            "n_points": int(len(aligned_stage)),
        }

        print(f"\nAligned {len(aligned_stage)} timesteps")
        print(f"  Period: {aligned_stage['datetime'].min()} to {aligned_stage['datetime'].max()}")
        print(f"  USGS gage datum: {gage_datum_ft:.2f} ft")
        print(f"  Modeled WSE range: {aligned_stage['modeled'].min():.2f} to {aligned_stage['modeled'].max():.2f} ft")
        print(f"  Observed WSE range: {aligned_stage['observed'].min():.2f} to {aligned_stage['observed'].max():.2f} ft")

        print("\n" + "="*60)
        print("STAGE / WSE VALIDATION METRICS")
        print("="*60)
        print(f"Bias: {stage_metrics['bias_ft']:.2f} ft")
        print(f"RMSE: {stage_metrics['rmse_ft']:.2f} ft")
        print(f"Mean Absolute Error: {stage_metrics['mae_ft']:.2f} ft")
        print(f"Peak Stage Error: {stage_metrics['peak_error_ft']:.2f} ft")
        print(f"Peak Timing Error: {stage_metrics['peak_time_error_hr']:.2f} hours")
else:
    print("Cannot calculate metrics - missing modeled or observed flow/stage data")
Text Only
Aligning modeled WSE and observed USGS stage timeseries...

Aligned 691 timesteps
  Period: 2020-12-22 05:00:00 to 2020-12-27 00:00:00
  USGS gage datum: 559.25 ft
  Modeled WSE range: 562.58 to 573.83 ft
  Observed WSE range: 565.92 to 571.03 ft

============================================================
STAGE / WSE VALIDATION METRICS
============================================================
Bias: 1.29 ft
RMSE: 2.67 ft
Mean Absolute Error: 2.01 ft
Peak Stage Error: 2.80 ft
Peak Timing Error: 37.67 hours
Python
# =============================================================================
# PLOT TIMESERIES COMPARISON
# =============================================================================
if aligned_flow is not None:
    fig = plot_timeseries_comparison(
        aligned_data=aligned_flow,
        metrics=flow_metrics,
        title=f"Flow Validation: {validation_meta.get('station_name', VALIDATION_GAUGE)}",
    )

    plt.tight_layout()
    plt.savefig(project_folder / "flow_validation.png", dpi=150)
    plt.show()

    print(f"\nFigure saved to: {project_folder / 'flow_validation.png'}")
elif aligned_stage is not None:
    fig, ax = plt.subplots(figsize=(12, 6))
    ax.plot(aligned_stage["datetime"], aligned_stage["observed"], "b-", linewidth=1.4, label="Observed USGS WSE")
    ax.plot(aligned_stage["datetime"], aligned_stage["modeled"], "r--", linewidth=1.4, label="Modeled nearest-cell WSE")
    ax.set_xlabel("Date")
    ax.set_ylabel("Water Surface Elevation (ft)")
    ax.set_title(f"Stage/WSE Validation: {validation_meta.get('station_name', VALIDATION_GAUGE)}")
    ax.grid(True, alpha=0.3)
    ax.legend(loc="best")

    metrics_text = (
        f"Datum = {stage_metrics['gage_datum_ft']:.2f} ft\n"
        f"Bias = {stage_metrics['bias_ft']:.2f} ft\n"
        f"RMSE = {stage_metrics['rmse_ft']:.2f} ft\n"
        f"Peak error = {stage_metrics['peak_error_ft']:.2f} ft"
    )
    ax.text(
        0.02, 0.98, metrics_text,
        transform=ax.transAxes,
        va="top",
        fontfamily="monospace",
        fontsize=9,
        bbox={"boxstyle": "round", "facecolor": "white", "alpha": 0.85},
    )

    plt.tight_layout()
    plt.savefig(project_folder / "stage_validation.png", dpi=150)
    plt.show()

    print(f"\nFigure saved to: {project_folder / 'stage_validation.png'}")
else:
    print("Cannot plot comparison - missing aligned data")

png

Text Only
Figure saved to: C:\Users\bill\.config\superpowers\worktrees\ras-commander-pr-docs\codex-pr251-timeseries-alignment\examples\example_projects\BaldEagleCrkMulti2D_914_historical\stage_validation.png
Python
# =============================================================================
# PLOT SCATTER / RESIDUAL COMPARISON
# =============================================================================
if aligned_flow is not None:
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    ax1, ax2 = axes

    ax1.scatter(aligned_flow["observed"], aligned_flow["modeled"], alpha=0.5, s=10)
    max_val = max(aligned_flow["observed"].max(), aligned_flow["modeled"].max())
    ax1.plot([0, max_val], [0, max_val], "r--", label="1:1 Line")
    ax1.set_xlabel("Observed Flow (cfs)")
    ax1.set_ylabel("Modeled Flow (cfs)")
    ax1.set_title(f"Flow Scatter Plot\nNSE={flow_metrics['nse']:.3f}, KGE={flow_metrics['kge']:.3f}")
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.set_aspect("equal")

    residuals = aligned_flow["modeled"] - aligned_flow["observed"]
    ax2.hist(residuals, bins=50, edgecolor="black", alpha=0.7)
    ax2.axvline(x=0, color="r", linestyle="--", label="Zero Error")
    ax2.axvline(x=residuals.mean(), color="g", linestyle="-", label=f"Mean: {residuals.mean():.0f} cfs")
    ax2.set_xlabel("Residual (Modeled - Observed) [cfs]")
    ax2.set_ylabel("Count")
    ax2.set_title("Residuals Distribution")
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig(project_folder / "flow_scatter.png", dpi=150)
    plt.show()

    print(f"\nFigure saved to: {project_folder / 'flow_scatter.png'}")
elif aligned_stage is not None:
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    ax1, ax2 = axes

    min_val = min(aligned_stage["observed"].min(), aligned_stage["modeled"].min())
    max_val = max(aligned_stage["observed"].max(), aligned_stage["modeled"].max())
    pad = max((max_val - min_val) * 0.05, 0.5)

    ax1.scatter(aligned_stage["observed"], aligned_stage["modeled"], alpha=0.5, s=10, color="darkorange")
    ax1.plot([min_val - pad, max_val + pad], [min_val - pad, max_val + pad], "k--", label="1:1 Line")
    ax1.set_xlabel("Observed WSE (ft)")
    ax1.set_ylabel("Modeled WSE (ft)")
    ax1.set_title(f"Stage/WSE Scatter\nRMSE={stage_metrics['rmse_ft']:.2f} ft, Bias={stage_metrics['bias_ft']:.2f} ft")
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.set_aspect("equal", adjustable="box")
    ax1.set_xlim(min_val - pad, max_val + pad)
    ax1.set_ylim(min_val - pad, max_val + pad)

    residuals = aligned_stage["modeled"] - aligned_stage["observed"]
    ax2.hist(residuals, bins=40, edgecolor="black", alpha=0.75, color="steelblue")
    ax2.axvline(x=0, color="r", linestyle="--", label="Zero Error")
    ax2.axvline(x=residuals.mean(), color="g", linestyle="-", label=f"Mean: {residuals.mean():.2f} ft")
    ax2.set_xlabel("Residual (Modeled - Observed) [ft]")
    ax2.set_ylabel("Count")
    ax2.set_title("Stage/WSE Residuals")
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig(project_folder / "stage_scatter.png", dpi=150)
    plt.show()

    print(f"\nFigure saved to: {project_folder / 'stage_scatter.png'}")
else:
    print("Cannot plot scatter - missing aligned data")

png

Text Only
Figure saved to: C:\Users\bill\.config\superpowers\worktrees\ras-commander-pr-docs\codex-pr251-timeseries-alignment\examples\example_projects\BaldEagleCrkMulti2D_914_historical\stage_scatter.png

Results Summary

Flow vs Stage Validation

This example first attempts flow validation from 1D cross-section output or 2D reference-line output. The active plan is a 2D-focused example and does not include a flow cutline/reference line at the downstream validation gauge, so the automated comparison falls back to stage validation:

  • Modeled value: 2D mesh-cell Water Surface Elevation at the cell nearest USGS 01548005.
  • Observed value: USGS gage height plus the reported USGS gage datum (gage_datum_ft).
  • Flow validation remains a model-development follow-up that requires adding a reviewed 2D reference line/cutline or computing flow from face velocities across a reviewed section.

Stage/WSE Metrics

Bias: mean modeled WSE minus observed WSE. Positive values indicate the model is high on average.

RMSE: root mean square stage error in feet.

Peak Stage Error: modeled peak WSE minus observed peak WSE.

Peak Timing Error: modeled peak time minus observed peak time in hours.

Expected Limitations

  1. Drainage Area Mismatch: The 2D model may not cover the entire HUC12 watershed, leading to missing runoff contributions from unmodeled areas.

  2. USGS Gauge Location: The validation gauge may be downstream of significant unmodeled tributaries.

  3. Precipitation Spatial Variability: AORC data at 4km resolution may not capture localized intense precipitation.

  4. Boundary Condition Uncertainty: USGS observations are used for validation context, while the active model run preserves the template model boundary hydrographs unless a later calibration step explicitly replaces them.

  5. Operations Data Constraint: Gate/control operations are preserved from the example model. Observed operations data are not included with this notebook.

  6. Point-Cell Stage Approximation: Nearest-cell WSE is appropriate for a demonstration fallback, but a professional validation should verify the gauge location, local datum, mesh resolution, and result extraction point.

Recommendations

  1. Add a Reviewed Flow Cutline: Add a 2D reference line at the validation gauge if flow metrics are needed.

  2. Expand Model Domain: Consider extending the 2D mesh to cover more of the watershed.

  3. Additional Validation Points: Add more USGS gauges if available within the model domain.

  4. Sensitivity Analysis: Test sensitivity to Manning's n, precipitation, and initial conditions.

  5. Multiple Events: Validate against multiple storm events to assess model robustness.

Python
# =============================================================================
# SUMMARY REPORT
# =============================================================================
print("="*70)
print("HISTORICAL EVENT VALIDATION SUMMARY")
print("="*70)
print(f"\nProject: {PROJECT_NAME}")
print(f"Storm Event: December 24-25, 2020")
print(f"Simulation Period: {SIM_START} to {SIM_END}")

print(f"\nData Sources:")
print(f"  Precipitation: AORC gridded data (4km, hourly)")
print(f"  Upstream BC: USGS {UPSTREAM_GAUGE}")
print(f"  Validation: USGS {VALIDATION_GAUGE}")

if HUC12_AVAILABLE and "coverage_pct" in dir():
    print(f"\nDrainage Coverage:")
    print(f"  Model covers {coverage_pct:.1f}% of HUC12 watershed")
    print(f"  Unmodeled area: {unmodeled_area_sqkm:.1f} sq km")

if flow_metrics is not None:
    print(f"\nFlow Validation Metrics:")
    print(f"  NSE: {flow_metrics['nse']:.3f}")
    print(f"  KGE: {flow_metrics['kge']:.3f}")
    print(f"  Peak Error: {flow_metrics['peak_error_pct']:.1f}%")
    print(f"  Volume Error: {flow_metrics['vol_error_pct']:.1f}%")
elif stage_metrics is not None:
    print(f"\nStage/WSE Validation Metrics:")
    print(f"  Gage datum: {stage_metrics['gage_datum_ft']:.2f} ft")
    print(f"  Bias: {stage_metrics['bias_ft']:.2f} ft")
    print(f"  RMSE: {stage_metrics['rmse_ft']:.2f} ft")
    print(f"  Mean Absolute Error: {stage_metrics['mae_ft']:.2f} ft")
    print(f"  Peak Stage Error: {stage_metrics['peak_error_ft']:.2f} ft")
    print(f"  Peak Timing Error: {stage_metrics['peak_time_error_hr']:.2f} hours")
    if stage_extraction_meta:
        print(
            f"  Extraction point: {stage_extraction_meta['mesh_name']} cell "
            f"{stage_extraction_meta['cell_id']} "
            f"({stage_extraction_meta['distance_ft']:.0f} ft from gauge)"
        )
    print("  Flow metrics not calculated because no reviewed flow cutline/reference line is present.")

print(f"\nOutput Files:")
if (project_folder / "coverage_analysis.png").exists():
    print(f"  - {project_folder / 'coverage_analysis.png'}")
if (project_folder / "flow_validation.png").exists():
    print(f"  - {project_folder / 'flow_validation.png'}")
if (project_folder / "flow_scatter.png").exists():
    print(f"  - {project_folder / 'flow_scatter.png'}")
if (project_folder / "stage_validation.png").exists():
    print(f"  - {project_folder / 'stage_validation.png'}")
if (project_folder / "stage_scatter.png").exists():
    print(f"  - {project_folder / 'stage_scatter.png'}")
if aorc_file and aorc_file.exists():
    print(f"  - {aorc_file}")

print("\n" + "="*70)
Text Only
======================================================================
HISTORICAL EVENT VALIDATION SUMMARY
======================================================================

Project: BaldEagleCrkMulti2D
Storm Event: December 24-25, 2020
Simulation Period: 2020-12-22 to 2020-12-27

Data Sources:
  Precipitation: AORC gridded data (4km, hourly)
  Upstream BC: USGS 01547500
  Validation: USGS 01548005

Stage/WSE Validation Metrics:
  Gage datum: 559.25 ft
  Bias: 1.29 ft
  RMSE: 2.67 ft
  Mean Absolute Error: 2.01 ft
  Peak Stage Error: 2.80 ft
  Peak Timing Error: 37.67 hours
  Extraction point: BaldEagleCr cell 10562 (97 ft from gauge)
  Flow metrics not calculated because no reviewed flow cutline/reference line is present.

Output Files:
  - C:\Users\bill\.config\superpowers\worktrees\ras-commander-pr-docs\codex-pr251-timeseries-alignment\examples\example_projects\BaldEagleCrkMulti2D_914_historical\stage_validation.png
  - C:\Users\bill\.config\superpowers\worktrees\ras-commander-pr-docs\codex-pr251-timeseries-alignment\examples\example_projects\BaldEagleCrkMulti2D_914_historical\stage_scatter.png
  - C:\Users\bill\.config\superpowers\worktrees\ras-commander-pr-docs\codex-pr251-timeseries-alignment\examples\example_projects\BaldEagleCrkMulti2D_914_historical\Precipitation\storm_20201224.nc

======================================================================

Next Steps for Modeler

Current State of Analysis

This notebook now completes an automated, event-aligned validation setup:

✅ Completed: - HUC12 drainage coverage analysis when optional watershed dependencies are available - AORC precipitation downloaded for Storm 12 (Dec 24-25, 2020) - USGS flow/stage records retrieved for validation context - Template plan and unsteady file cloned for the 2020 event - Cloned plan simulation window aligned to the AORC event window - Cloned unsteady file configured to use Precipitation/storm_20201224.nc - Existing model boundary and gate/control time series preserved from the template model - Nearest-cell 2D WSE extracted and compared to observed USGS stage converted with the reported gage datum

⚠️ Constraints and Assumptions: 1. USGS observations are not injected into the active unsteady file in this notebook; they are used for comparison and validation context. 2. Existing upstream/downstream boundary conditions remain template-model assumptions unless a later calibration step replaces them. 3. Gate/control operations are preserved from the example project because observed operations are not included. 4. Additional tributary inflows may be needed for a professional validation model if the downstream gauge includes drainage area not represented by the active mesh. 5. Flow metrics are not calculated until a reviewed 2D reference line/cutline or face-flow integration method is added at the validation location. 6. The stage comparison uses the nearest mesh-cell WSE as a demonstration fallback; a professional validation should review the extraction location, datum, and mesh resolution.

Phase 1: Boundary-Condition Calibration Follow-Up

To convert this from event-aligned forcing validation into a calibrated historical event model, update the active cloned unsteady file with observed hydrographs where appropriate:

  1. Replace the existing upstream flow hydrograph with USGS 01547500 flow data if that gauge is accepted as the upstream boundary condition.
  2. Evaluate whether the downstream normal-depth boundary should remain a hydraulic boundary or be replaced with a stage hydrograph from an appropriate gauge.
  3. Add or estimate lateral inflows for material tributaries that enter between the upstream and validation gauges.
  4. Add a reviewed 2D reference line or cutline at the validation gauge if flow metrics are needed.
  5. Re-run the model and compare modeled results against the validation gauge after boundary-condition updates.

Phase 2: Enhanced Validation with Lateral Inflows (Future)

Add new boundary conditions for lateral tributaries to improve validation.

2.1 Required New BCs

Tributary Gauge DA (sq mi) Data Available Priority
Spring Creek 01547100 142 Flow (IV) ✓ HIGH
Marsh Creek 01547700 44 Flow (IV) ✓ HIGH
Beech Creek 01547980 170 No data Low
Fishing Creek 01548079 180 No data Low

Combined gaged inflow: 339 + 142 + 44 = 525 sq mi (vs 562 sq mi at downstream gauge) Ungauged gap: 37 sq mi (6.6% of downstream drainage)

2.2 Geometry File Edits Required

For each new BC, you need to:

  1. Create SA/2D Area Conn in geometry file:
  2. Add connection line between external boundary and 2D mesh
  3. Define as "SA/2D Area Conn" type

  4. Specify connection cells

  5. Define External SA (storage area):

  6. Create storage area for lateral inflow

  7. Connect to 2D mesh via SA/2D Area Conn

  8. Add BC Reference in unsteady file:

  9. Reference the new SA connection
  10. Add flow hydrograph table

Current Limitation: This requires manual geometry editing or HEC-RAS GUI operations.

2.3 Future Automation Opportunity

There is not currently a public ras-commander API that creates new lateral-inflow boundary geometry directly from a gauge coordinate. For this notebook, treat the tributary inflow boundaries as a manual geometry-authoring task that should be reviewed in HEC-RAS/RAS Mapper before simulation.

A future helper could support this workflow by:

  1. identifying candidate 2D mesh faces near a gauge location;
  2. building a defensible connection alignment outside the mesh;
  3. trimming or smoothing the alignment to avoid mesh corners;
  4. adding the SA/2D Area Connection or other required geometry records;
  5. adding the corresponding unsteady-flow boundary reference; and
  6. validating the edited geometry before use.

Until that capability exists, do not assume a one-call gauge-to-boundary API is available. Use the existing geometry and unsteady-flow editing APIs only where they explicitly support the required HEC-RAS record type, and perform geometry edits with normal modeler review.

Phase 3: Multi-Gauge Validation Network

Once lateral BCs are added, validate at multiple points:

Validation Point Gauge DA at Gauge DA Upstream BCs Coverage
Downstream 01548005 562 sq mi 525 sq mi 93.4%
Midstream 01548000 559 sq mi 525 sq mi 93.9%
Spring Creek 01547100 142 sq mi 142 sq mi 100% (at BC)

Validation Approach: - Extract results at reviewed cross sections, reference lines, or mesh cells near each gauge - Compare modeled vs observed stage and/or flow depending on available outputs - Calculate metrics for each validation point - Assess spatial performance (upstream vs downstream)

Phase 4: Multi-Event Calibration (Future)

Use AORC catalog to test multiple storms:

From the 2020 AORC catalog: 1. Storm 12 (Dec 24-25): 2.72 in - largest event 2. Storm 5 (Apr 30-May 1): 2.12 in - spring conditions 3. Storm 10 (Nov 11): 1.99 in - high intensity 4. Storm 9 (Oct 29-30): 1.74 in - fall conditions

Multi-event validation benefits: - Test across different magnitudes - Assess parameter transferability - Identify seasonal biases - Build confidence in calibration

Tools and References

Automation Functions Available: - RasUnsteady.set_gridded_precipitation() - Configure AORC NetCDF - retrieve_flow_data() / retrieve_stage_data() - Retrieve USGS observations - calculate_all_metrics() - Comprehensive flow validation metrics when modeled flow is available - plot_timeseries_comparison() - Flow comparison plots when modeled flow is available

Documentation: - BC Configuration: .claude/outputs/general-purpose/2025-12-30-bc-configuration-workflow.md - Gauge Analysis: .claude/outputs/general-purpose/2025-12-29-gauge-data-availability.md - AORC Workflow: .claude/outputs/general-purpose/2025-12-29-aorc-workflow-research.md - Complete Summary: .claude/outputs/general-purpose/2025-12-30-validation-workflow-summary.md

Estimated Time to Complete: - Phase 1 (Minimum viable): 2-3 hours - Phase 2 (Add lateral BCs): +3-5 hours (pending automation feature) - Phase 3 (Multi-gauge validation): +2-3 hours - Phase 4 (Multi-event): +4-6 hours per additional storm

Recommended Approach: Start with Phase 1, validate the workflow works, then expand.

CLB Engineering Corporation  ·  LLM Forward Engineering
RAS Commander Documentation is a free and open-source project maintained by CLB Engineering Corporation. For agencies and firms modernizing H&H workflows with LLM Forward approaches, contact CLB to partner with the engineers who wrote the automation.