Skip to content

SA/2D Connection Authoring

Create, inspect, delete, and re-create SA/2D area connections in a real HEC-RAS geometry file. This notebook demonstrates an engineering-oriented workflow for the Dam connection in BaldEagleDamBrk.g13:

  • inspect the connection line, crest profile, gates, terrain, mesh, and boundary context
  • sample the project terrain raster along the connection polyline
  • re-create the connection and write a terrain-derived weir crest profile
  • raise the minimum terrain-derived crest by 0.5 ft and floor lower points to that value
  • verify the result with before/after profile and map views

Primary APIs used:

  • GeomLateral.get_connections() - list SA/2D connections
  • GeomLateral.get_connection_line_coords() - read connection polyline coordinates
  • GeomLateral.get_connection_profile() - read weir crest station/elevation profile
  • GeomLateral.get_connection_gates() - read gate definitions
  • GeomLateral.delete_connection() - remove a connection block
  • GeomLateral.set_connection() - create or replace a connection block
  • GeomLateral.set_connection_profile() - write the sampled/modified weir crest profile
  • GeomLateral.set_connection_gates() - write gate definitions
  • HdfMesh and HdfBndry - read mesh cells, mesh perimeters, and BC lines from the geometry HDF

The terrain profile in this notebook is sampled directly from the extracted example project's Terrain/Terrain50.dtm_20ft.tif raster. RasExamples.extract_project() provides that terrain raster, the .rasmap, and the compiled geometry HDF used for the map context.

Development Mode

Python
#!pip install --upgrade ras-commander
Python
# =============================================================================
# DEVELOPMENT MODE TOGGLE
# =============================================================================
USE_LOCAL_SOURCE = False

if USE_LOCAL_SOURCE:
    import sys
    from pathlib import Path

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

from pathlib import Path
import os
import logging
import warnings

import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from matplotlib.colors import Normalize
from matplotlib.lines import Line2D
from matplotlib.patches import Rectangle
import numpy as np
import pandas as pd
import rasterio
from rasterio.enums import Resampling
from rasterio.windows import bounds as window_bounds
from rasterio.windows import from_bounds
from shapely.geometry import LineString, box

from IPython.display import display

from ras_commander import RasExamples
from ras_commander.geom import GeomLateral
from ras_commander.hdf import HdfBndry, HdfMesh

warnings.filterwarnings("ignore", category=FutureWarning)
logging.getLogger("ras_commander").setLevel(logging.CRITICAL)

pd.set_option("display.max_columns", None)
pd.set_option("display.max_colwidth", 120)

plt.rcParams.update({
    "figure.dpi": 120,
    "axes.grid": True,
    "grid.alpha": 0.22,
    "axes.spines.top": False,
    "axes.spines.right": False,
})


def connection_line(coords_df):
    return LineString(coords_df[["X", "Y"]].to_numpy())


def expanded_bounds(bounds, fraction=0.25, min_pad=100.0):
    minx, miny, maxx, maxy = bounds
    width = max(maxx - minx, min_pad)
    height = max(maxy - miny, min_pad)
    pad_x = max(width * fraction, min_pad)
    pad_y = max(height * fraction, min_pad)
    return (minx - pad_x, miny - pad_y, maxx + pad_x, maxy + pad_y)


def sample_points_along_line(coords_df, spacing):
    line = connection_line(coords_df)
    stations = np.arange(0.0, line.length, spacing)
    if len(stations) == 0 or not np.isclose(stations[-1], line.length):
        stations = np.r_[stations, line.length]
    points = [line.interpolate(float(sta)) for sta in stations]
    return pd.DataFrame({
        "Station": np.round(stations, 3),
        "X": [pt.x for pt in points],
        "Y": [pt.y for pt in points],
    })


def sample_terrain_profile(coords_df, terrain_raster, spacing):
    samples = sample_points_along_line(coords_df, spacing)
    xy = list(zip(samples["X"], samples["Y"]))
    elevations = []
    with rasterio.open(terrain_raster) as src:
        for value in src.sample(xy, masked=True):
            elev = value[0]
            elevations.append(np.nan if np.ma.is_masked(elev) else float(elev))

    samples["Elevation"] = elevations
    if samples["Elevation"].isna().any():
        samples["Elevation"] = samples["Elevation"].interpolate(limit_direction="both")
    if samples["Elevation"].isna().any():
        raise ValueError("Terrain sampling produced unresolved NoData elevations.")
    return samples[["Station", "Elevation", "X", "Y"]]


def hillshade(elev_ma, azimuth=315, altitude=45):
    arr = np.ma.filled(elev_ma.astype(float), np.nan)
    if np.isnan(arr).all():
        return np.ma.masked_all(arr.shape)
    arr = np.where(np.isnan(arr), np.nanmean(arr), arr)
    dy, dx = np.gradient(arr)
    slope = np.pi / 2.0 - np.arctan(np.hypot(dx, dy))
    aspect = np.arctan2(-dx, dy)
    az = np.deg2rad(azimuth)
    alt = np.deg2rad(altitude)
    shaded = (
        np.sin(alt) * np.sin(slope)
        + np.cos(alt) * np.cos(slope) * np.cos(az - aspect)
    )
    return np.clip((shaded + 1.0) / 2.0, 0.0, 1.0)


def read_terrain_window(terrain_raster, bounds, max_pixels=900):
    left, bottom, right, top = bounds
    with rasterio.open(terrain_raster) as src:
        window = from_bounds(left, bottom, right, top, transform=src.transform)
        window = window.round_offsets().round_lengths()
        scale = max(window.width / max_pixels, window.height / max_pixels, 1.0)
        out_height = max(1, int(window.height / scale))
        out_width = max(1, int(window.width / scale))
        elev = src.read(
            1,
            window=window,
            out_shape=(out_height, out_width),
            masked=True,
            boundless=True,
            resampling=Resampling.bilinear,
        )
        left, bottom, right, top = window_bounds(window, src.transform)
        extent = (left, right, bottom, top)
    return elev, extent


def add_north_arrow(ax):
    ax.annotate(
        "N",
        xy=(0.94, 0.90),
        xytext=(0.94, 0.78),
        xycoords="axes fraction",
        ha="center",
        va="center",
        fontsize=10,
        arrowprops=dict(arrowstyle="-|>", color="black", linewidth=1.5),
        bbox=dict(boxstyle="round,pad=0.18", facecolor="white", edgecolor="0.7", alpha=0.8),
        zorder=20,
    )


def nice_scale_length(raw_length):
    if raw_length <= 0:
        return 100.0
    exponent = np.floor(np.log10(raw_length))
    base = raw_length / (10 ** exponent)
    step = 1 if base < 1.5 else 2 if base < 3.5 else 5 if base < 7.5 else 10
    return float(step * (10 ** exponent))


def add_scale_bar(ax):
    xmin, xmax = ax.get_xlim()
    ymin, ymax = ax.get_ylim()
    width = xmax - xmin
    height = ymax - ymin
    length = nice_scale_length(width * 0.18)
    x0 = xmin + width * 0.06
    y0 = ymin + height * 0.07
    ax.plot([x0, x0 + length], [y0, y0], color="black", linewidth=3, zorder=20)
    ax.plot([x0, x0], [y0 - height * 0.01, y0 + height * 0.01], color="black", linewidth=1.5, zorder=20)
    ax.plot([x0 + length, x0 + length], [y0 - height * 0.01, y0 + height * 0.01], color="black", linewidth=1.5, zorder=20)
    label = f"{length:,.0f} ft" if length < 5280 else f"{length / 5280:.1f} mi"
    ax.text(
        x0 + length / 2,
        y0 + height * 0.025,
        label,
        ha="center",
        va="bottom",
        fontsize=8,
        bbox=dict(facecolor="white", edgecolor="0.7", alpha=0.8, pad=2),
        zorder=20,
    )


def dedupe_legend(ax, loc="upper right", fontsize=8):
    handles, labels = ax.get_legend_handles_labels()
    seen = {}
    for handle, label in zip(handles, labels):
        if label and label not in seen:
            seen[label] = handle
    if seen:
        ax.legend(seen.values(), seen.keys(), loc=loc, fontsize=fontsize, frameon=True)


def label_lines(ax, gdf, column, bounds, fontsize=8):
    if gdf.empty:
        return
    view = box(*bounds)
    for _, row in gdf.iterrows():
        geom = row.geometry
        if geom is None or geom.is_empty or not geom.intersects(view):
            continue
        label_pt = geom.interpolate(0.5, normalized=True)
        ax.text(
            label_pt.x,
            label_pt.y,
            str(row[column]),
            fontsize=fontsize,
            ha="center",
            va="center",
            bbox=dict(facecolor="white", edgecolor="0.7", alpha=0.75, pad=1.5),
            zorder=15,
        )


def plot_context_map(
    ax,
    terrain_raster,
    mesh_areas,
    mesh_cells,
    bc_lines,
    bounds,
    title,
    show_cells=True,
):
    elev, extent = read_terrain_window(terrain_raster, bounds)
    valid = elev.compressed()
    if len(valid):
        vmin, vmax = np.percentile(valid, [2, 98])
    else:
        vmin, vmax = None, None
    ax.imshow(
        elev,
        extent=extent,
        origin="upper",
        cmap="gist_earth",
        vmin=vmin,
        vmax=vmax,
        alpha=0.82,
        zorder=0,
    )
    ax.imshow(hillshade(elev), extent=extent, origin="upper", cmap="gray", alpha=0.22, zorder=1)
    if len(valid):
        left_e, right_e, bottom_e, top_e = extent
        xs = np.linspace(left_e, right_e, elev.shape[1])
        ys = np.linspace(top_e, bottom_e, elev.shape[0])
        levels = np.linspace(vmin, vmax, 9)
        ax.contour(
            xs,
            ys,
            np.ma.filled(elev, np.nan),
            levels=levels,
            colors="white",
            linewidths=0.3,
            alpha=0.45,
            zorder=2,
        )

    left, bottom, right, top = bounds
    if show_cells and not mesh_cells.empty:
        cells_view = mesh_cells.cx[left:right, bottom:top]
        if not cells_view.empty:
            cells_view.boundary.plot(ax=ax, color="#4a5568", linewidth=0.28, alpha=0.50, zorder=3, label="2D mesh cells")

    if not mesh_areas.empty:
        mesh_areas.boundary.plot(ax=ax, color="#1f4e79", linewidth=1.2, alpha=0.95, zorder=4, label="2D flow area")
        for _, row in mesh_areas.iterrows():
            pt = row.geometry.representative_point()
            if left <= pt.x <= right and bottom <= pt.y <= top:
                ax.text(
                    pt.x,
                    pt.y,
                    row["mesh_name"],
                    fontsize=9,
                    weight="bold",
                    ha="center",
                    va="center",
                    bbox=dict(facecolor="white", edgecolor="#1f4e79", alpha=0.7, pad=2),
                    zorder=14,
                )

    if not bc_lines.empty:
        bc_view = bc_lines.cx[left:right, bottom:top]
        if not bc_view.empty:
            bc_view.plot(ax=ax, color="#2b6cb0", linewidth=2.0, linestyle="--", zorder=6, label="BC line")
            label_lines(ax, bc_view, "Name", bounds)

    ax.set_xlim(left, right)
    ax.set_ylim(bottom, top)
    ax.set_aspect("equal", adjustable="box")
    ax.set_title(title)
    ax.set_xlabel("Easting (ft)")
    ax.set_ylabel("Northing (ft)")
    add_north_arrow(ax)
    add_scale_bar(ax)


def add_gate_rectangles(ax, gates, facecolor="#e17c05", edgecolor="#9a3412", label="Gate opening"):
    label_used = False
    if gates.empty:
        return
    for _, gate in gates.iterrows():
        stations = gate.get("OpeningStations", []) or []
        width = float(gate.get("Width", 0.0) or 0.0)
        height = float(gate.get("Height", 0.0) or 0.0)
        invert = float(gate.get("InvertElevation", np.nan))
        if width <= 0 or height <= 0 or np.isnan(invert):
            continue
        for sta in stations:
            rect = Rectangle(
                (float(sta) - width / 2.0, invert),
                width,
                height,
                facecolor=facecolor,
                edgecolor=edgecolor,
                linewidth=1.1,
                alpha=0.28,
                label=label if not label_used else None,
                zorder=8,
            )
            ax.add_patch(rect)
            ax.text(
                float(sta),
                invert + height,
                gate["GateName"],
                fontsize=7,
                ha="center",
                va="bottom",
                color=edgecolor,
                rotation=90,
                zorder=9,
            )
            label_used = True


def set_profile_ylim(ax, profiles, gates=None):
    if gates is None:
        gates = pd.DataFrame()
    values = []
    for profile in profiles:
        values.extend(profile["Elevation"].dropna().tolist())
    if not gates.empty:
        for _, gate in gates.iterrows():
            invert = float(gate.get("InvertElevation", np.nan))
            height = float(gate.get("Height", 0.0) or 0.0)
            if not np.isnan(invert):
                values.extend([invert, invert + height])
    if values:
        ymin, ymax = min(values), max(values)
        pad = max((ymax - ymin) * 0.08, 2.0)
        ax.set_ylim(ymin - pad, ymax + pad)


def add_gate_station_markers(ax, gates):
    if gates.empty:
        return
    top = ax.get_ylim()[1]
    for _, gate in gates.iterrows():
        stations = gate.get("OpeningStations", []) or []
        for sta in stations:
            ax.axvline(float(sta), color="#9a3412", linewidth=0.9, alpha=0.45, zorder=3)
            ax.text(
                float(sta),
                top,
                gate["GateName"],
                fontsize=7,
                ha="center",
                va="top",
                color="#9a3412",
                rotation=90,
                zorder=9,
            )


def plot_gate_band(ax, gates):
    if gates.empty:
        ax.text(0.5, 0.5, "No gates", transform=ax.transAxes, ha="center", va="center")
        return
    add_gate_rectangles(ax, gates)
    values = []
    for _, gate in gates.iterrows():
        invert = float(gate.get("InvertElevation", np.nan))
        height = float(gate.get("Height", 0.0) or 0.0)
        if not np.isnan(invert):
            values.extend([invert, invert + height])
    ymin, ymax = min(values), max(values)
    pad = max((ymax - ymin) * 0.25, 2.0)
    ax.set_ylim(ymin - pad, ymax + pad)
    ax.set_ylabel("Gate elev. (ft)")
    ax.set_xlabel("Station along connection (ft)")
    ax.set_title("Gate Opening Elevations")
    dedupe_legend(ax, loc="upper right", fontsize=7)


def plot_profile(ax, terrain_profile, crest_profile, gates, title, crest_label, show_gate_rectangles=True):
    ax.plot(
        terrain_profile["Station"],
        terrain_profile["Elevation"],
        color="0.55",
        linewidth=1.4,
        linestyle="--",
        label="Terrain profile",
        zorder=2,
    )
    ax.plot(
        crest_profile["Station"],
        crest_profile["Elevation"],
        color="#1f4e79",
        linewidth=2.1,
        label=crest_label,
        zorder=4,
    )
    if show_gate_rectangles:
        add_gate_rectangles(ax, gates)
    ax.set_title(title)
    ax.set_xlabel("Station along connection (ft)")
    ax.set_ylabel("Elevation (ft)")
    set_profile_ylim(ax, [terrain_profile, crest_profile], gates if show_gate_rectangles else None)
    dedupe_legend(ax, loc="lower right")


def save_figure(fig, filename):
    FIGURE_DIR.mkdir(parents=True, exist_ok=True)
    output_path = FIGURE_DIR / filename
    fig.savefig(output_path, dpi=180)
    print(f"Saved figure: {output_path}")


import ras_commander

print(f"Loaded: {ras_commander.__file__}")
Text Only
PIP PACKAGE MODE: loading installed ras-commander


Loaded: <workspace>\ras_commander\__init__.py

Parameters

Python
PROJECT_NAME = "BaldEagleCrkMulti2D"
PROJECT_SUFFIX = "214_connection_authoring"
GEOM_FILE_NAME = "BaldEagleDamBrk.g13"
TARGET_CONNECTION = "Dam"

TERRAIN_RASTER_NAME = "Terrain50.dtm_20ft.tif"
TERRAIN_SAMPLE_SPACING_FT = 25.0
CREST_IMPROVEMENT_FT = 0.5

cwd = Path.cwd()
REPO_ROOT = cwd if (cwd / "ras_commander").exists() else cwd.parent
WORK_ROOT = REPO_ROOT / "working" / "connection_authoring"
WORK_ROOT.mkdir(parents=True, exist_ok=True)
FIGURE_DIR = Path(
    os.environ.get(
        "RAS_COMMANDER_NOTEBOOK_FIGURE_DIR",
        str(WORK_ROOT / "figures"),
    )
)
FIGURE_DIR.mkdir(parents=True, exist_ok=True)

print(f"Working folder: {WORK_ROOT}")
print(f"Figure folder: {FIGURE_DIR}")
Text Only
Working folder: <workspace>\working\connection_authoring
Figure folder: <workspace>\working\connection_authoring\figures

Extract And Inspect The Project

Python
project_path = RasExamples.extract_project(
    PROJECT_NAME,
    output_path=WORK_ROOT,
    suffix=PROJECT_SUFFIX,
)
geom_file = project_path / GEOM_FILE_NAME
geom_hdf = project_path / f"{GEOM_FILE_NAME}.hdf"
rasmap_file = project_path / "BaldEagleDamBrk.rasmap"
terrain_raster = project_path / "Terrain" / TERRAIN_RASTER_NAME

print(f"Project path: {project_path}")
print(f"Geometry file: {geom_file.name}")
print(f"Geometry HDF: {geom_hdf.name}")
print(f"RASMapper file: {rasmap_file.name}")
print(f"Terrain raster: {terrain_raster.relative_to(project_path)}")

assert geom_file.exists()
assert geom_hdf.exists()
assert rasmap_file.exists()
assert terrain_raster.exists()

connections = GeomLateral.get_connections(geom_file)
display_cols = ["Name", "Type", "From", "To", "NumPoints", "LinePoints", "HasGate"]
display(connections[display_cols])

mesh_areas = HdfMesh.get_mesh_areas(geom_hdf)
mesh_cells = HdfMesh.get_mesh_cell_polygons(geom_hdf)
bc_lines = HdfBndry.get_bc_lines(geom_hdf)

connection_lines = {}
for name in connections["Name"]:
    try:
        connection_lines[name] = GeomLateral.get_connection_line_coords(geom_file, name)
    except ValueError:
        pass

context_summary = pd.DataFrame([{
    "Connections": len(connections),
    "2D Flow Areas": len(mesh_areas),
    "2D Mesh Cells": len(mesh_cells),
    "BC Lines": len(bc_lines),
    "Terrain Source": str(terrain_raster.relative_to(project_path)),
    "Mesh Source": geom_hdf.name,
}])
display(context_summary)
Text Only
Project path: <workspace>\working\connection_authoring\BaldEagleCrkMulti2D_214_connection_authoring
Geometry file: BaldEagleDamBrk.g13
Geometry HDF: BaldEagleDamBrk.g13.hdf
RASMapper file: BaldEagleDamBrk.rasmap
Terrain raster: Terrain\Terrain50.dtm_20ft.tif
Name Type From To NumPoints LinePoints HasGate
0 Dam SA to 2D Reservoir Pool BaldEagleCr 6 18 True
1 Lower Levee 2D to 2D BaldEagleCr BaldEagleCr 94 89 False
2 Middle Levee 2D to 2D BaldEagleCr BaldEagleCr 130 122 False
3 Upper Levee 2D to 2D BaldEagleCr BaldEagleCr 101 99 False
Connections 2D Flow Areas 2D Mesh Cells BC Lines Terrain Source Mesh Source
0 4 1 14364 2 Terrain\Terrain50.dtm_20ft.tif BaldEagleDamBrk.g13.hdf

Model Context And Visualization Targets

Before modifying the connection, an engineer needs enough spatial and hydraulic context to interpret the profile:

  • target connection alignment and adjacent storage/2D flow areas
  • 2D mesh perimeter and local cells crossed by the connection
  • downstream boundary condition locations
  • underlying terrain along and around the connection
  • weir crest station/elevation profile and gate openings

The overview below locates the target connection in the model. The detailed inspection maps that follow zoom to the connection line with a 25 percent margin around its extents.

Python
target_coords = connection_lines[TARGET_CONNECTION]
target_line = connection_line(target_coords)
overview_bounds = expanded_bounds(mesh_areas.total_bounds, fraction=0.04, min_pad=1000.0)

fig, ax = plt.subplots(figsize=(11, 7))
plot_context_map(
    ax,
    terrain_raster,
    mesh_areas,
    mesh_cells,
    bc_lines,
    overview_bounds,
    "Model Overview: Mesh, BC Lines, And SA/2D Connections",
    show_cells=False,
)

for name, coords in connection_lines.items():
    color = "#b91c1c" if name == TARGET_CONNECTION else "#7c2d12"
    linewidth = 2.8 if name == TARGET_CONNECTION else 1.4
    label = f"Target connection: {name}" if name == TARGET_CONNECTION else "Other SA/2D connections"
    ax.plot(coords["X"], coords["Y"], color=color, linewidth=linewidth, alpha=0.95, zorder=8, label=label)
    mid = connection_line(coords).interpolate(0.5, normalized=True)
    ax.text(
        mid.x,
        mid.y,
        name,
        fontsize=8,
        ha="center",
        va="center",
        bbox=dict(facecolor="white", edgecolor=color, alpha=0.75, pad=1.5),
        zorder=16,
    )

dedupe_legend(ax, loc="lower left")
fig.tight_layout()
save_figure(fig, "214_connection_authoring_model_overview.png")
plt.show()
Text Only
Saved figure: <workspace>\working\connection_authoring\figures\214_connection_authoring_model_overview.png

jpeg

Inspect The Target Connection

Python
original_coords = GeomLateral.get_connection_line_coords(geom_file, TARGET_CONNECTION)
original_profile = GeomLateral.get_connection_profile(geom_file, TARGET_CONNECTION)
original_gates = GeomLateral.get_connection_gates(geom_file, TARGET_CONNECTION)
terrain_profile = sample_terrain_profile(
    original_coords,
    terrain_raster,
    TERRAIN_SAMPLE_SPACING_FT,
)

print(f"Connection line points: {len(original_coords)}")
print(f"Original weir crest points: {len(original_profile)}")
print(f"Terrain profile samples: {len(terrain_profile)} at {TERRAIN_SAMPLE_SPACING_FT:g} ft spacing")
print(f"Gates: {len(original_gates)}")

gate_cols = ["GateName", "Width", "Height", "InvertElevation", "GateCoefficient", "NumOpenings", "OpeningStations"]
profile_stats = pd.DataFrame([{
    "Source": "Existing weir crest",
    "Points": len(original_profile),
    "Min Elev": original_profile["Elevation"].min(),
    "Max Elev": original_profile["Elevation"].max(),
    "Min Station": original_profile.loc[original_profile["Elevation"].idxmin(), "Station"],
}, {
    "Source": "Terrain raster sample",
    "Points": len(terrain_profile),
    "Min Elev": terrain_profile["Elevation"].min(),
    "Max Elev": terrain_profile["Elevation"].max(),
    "Min Station": terrain_profile.loc[terrain_profile["Elevation"].idxmin(), "Station"],
}])
display(profile_stats.round(3))
display(original_gates[gate_cols])
Text Only
Connection line points: 18
Original weir crest points: 6
Terrain profile samples: 298 at 25 ft spacing
Gates: 1
Source Points Min Elev Max Elev Min Station
0 Existing weir crest 6 657.000 683.0 2250.0
1 Terrain raster sample 298 647.656 692.5 2275.0
GateName Width Height InvertElevation GateCoefficient NumOpenings OpeningStations
0 Gate #1 7.0 15.0 590.0 0.65 2 [5745.0, 5765.0]

Weir Crest Profile With Terrain And Gates

Python
fig, (ax, gate_ax) = plt.subplots(
    2,
    1,
    figsize=(12, 6.2),
    sharex=True,
    gridspec_kw={"height_ratios": [4.0, 1.15]},
)
plot_profile(
    ax,
    terrain_profile,
    original_profile,
    original_gates,
    f"{TARGET_CONNECTION}: Existing Weir Crest Against Terrain",
    "Existing weir crest",
    show_gate_rectangles=False,
)
add_gate_station_markers(ax, original_gates)
plot_gate_band(gate_ax, original_gates)
fig.tight_layout()
save_figure(fig, "214_connection_authoring_existing_profile.png")
plt.show()
Text Only
Saved figure: <workspace>\working\connection_authoring\figures\214_connection_authoring_existing_profile.png

jpeg

Connection Map View

Python
connection_bounds = expanded_bounds(target_line.bounds, fraction=0.25, min_pad=250.0)

fig, ax = plt.subplots(figsize=(11, 7))
plot_context_map(
    ax,
    terrain_raster,
    mesh_areas,
    mesh_cells,
    bc_lines,
    connection_bounds,
    f"{TARGET_CONNECTION}: Terrain, Mesh Cells, And Connection Alignment",
    show_cells=True,
)
ax.plot(original_coords["X"], original_coords["Y"], color="#b91c1c", linewidth=2.8, zorder=10, label="Connection line")
ax.scatter(original_coords["X"], original_coords["Y"], s=16, color="#fee2e2", edgecolor="#b91c1c", linewidth=0.6, zorder=11)

dedupe_legend(ax, loc="upper left")
fig.tight_layout()
save_figure(fig, "214_connection_authoring_connection_map.png")
plt.show()
Text Only
Saved figure: <workspace>\working\connection_authoring\figures\214_connection_authoring_connection_map.png

jpeg

Delete The Connection

Python
GeomLateral.delete_connection(geom_file, TARGET_CONNECTION, create_backup=True)

after_delete = GeomLateral.get_connections(geom_file)
display(after_delete[["Name", "Type", "From", "To"]])

assert TARGET_CONNECTION not in after_delete["Name"].values
assert len(after_delete) == len(connections) - 1
print(f"Deleted '{TARGET_CONNECTION}'. {len(after_delete)} connections remain.")
Name Type From To
0 Lower Levee 2D to 2D BaldEagleCr BaldEagleCr
1 Middle Levee 2D to 2D BaldEagleCr BaldEagleCr
2 Upper Levee 2D to 2D BaldEagleCr BaldEagleCr
Text Only
Deleted 'Dam'. 3 connections remain.

Re-Create The Connection With Modified Parameters

Python
coords_tuples = list(zip(
    original_coords["X"].tolist(),
    original_coords["Y"].tolist(),
))

dam_row = connections.loc[connections["Name"] == TARGET_CONNECTION].iloc[0]

GeomLateral.set_connection(
    geom_file,
    TARGET_CONNECTION,
    coords_tuples,
    upstream_area=dam_row["From"],
    downstream_area=dam_row["To"],
    routing_type=1,
    weir_width=120.0,
    weir_coef=2.6,
    overflow_method_2d=True,
    create_backup=False,
)

after_recreate = GeomLateral.get_connections(geom_file)
display(after_recreate[["Name", "Type", "From", "To", "NumPoints", "LinePoints"]])

assert TARGET_CONNECTION in after_recreate["Name"].values
assert len(after_recreate) == len(connections)
print(f"Re-created '{TARGET_CONNECTION}' with modified weir coefficient and width.")
Name Type From To NumPoints LinePoints
0 Lower Levee 2D to 2D BaldEagleCr BaldEagleCr 94 89
1 Middle Levee 2D to 2D BaldEagleCr BaldEagleCr 130 122
2 Upper Levee 2D to 2D BaldEagleCr BaldEagleCr 101 99
3 Dam SA to 2D Reservoir Pool BaldEagleCr 2 18
Text Only
Re-created 'Dam' with modified weir coefficient and width.

Derive And Modify The Weir Crest From Terrain

Python
terrain_derived_profile = terrain_profile[["Station", "Elevation", "X", "Y"]].copy()
terrain_min_idx = terrain_derived_profile["Elevation"].idxmin()
terrain_min_elev = float(terrain_derived_profile.loc[terrain_min_idx, "Elevation"])
new_min_elev = terrain_min_elev + CREST_IMPROVEMENT_FT

modified_profile = terrain_derived_profile.copy()
modified_profile["TerrainElevation"] = terrain_derived_profile["Elevation"]
modified_profile["Elevation"] = np.maximum(modified_profile["TerrainElevation"], new_min_elev)
modified_profile["Delta"] = modified_profile["Elevation"] - modified_profile["TerrainElevation"]

GeomLateral.set_connection_profile(
    geom_file,
    TARGET_CONNECTION,
    modified_profile[["Station", "Elevation"]],
    create_backup=False,
)

updated_profile = GeomLateral.get_connection_profile(geom_file, TARGET_CONNECTION)
raised_points = modified_profile[modified_profile["Delta"] > 0].copy()

display(pd.DataFrame([{
    "Terrain Min Station": modified_profile.loc[terrain_min_idx, "Station"],
    "Terrain Min Elev": terrain_min_elev,
    "New Minimum Crest Elev": new_min_elev,
    "Improvement": CREST_IMPROVEMENT_FT,
    "Points Raised": len(raised_points),
    "Profile Points Written": len(updated_profile),
}]).round(3))
display(raised_points[["Station", "TerrainElevation", "Elevation", "Delta"]].round(3).head(12))

assert len(updated_profile) == len(modified_profile)
assert np.isclose(updated_profile["Elevation"].min(), new_min_elev, atol=0.001)
assert (modified_profile["Elevation"] >= new_min_elev - 1e-9).all()
print(
    f"Terrain-derived profile written. Minimum crest raised from "
    f"{terrain_min_elev:.2f} ft to {new_min_elev:.2f} ft."
)
Terrain Min Station Terrain Min Elev New Minimum Crest Elev Improvement Points Raised Profile Points Written
0 2275.0 647.656 648.156 0.5 16 298
Station TerrainElevation Elevation Delta
91 2275.0 647.656 648.156 0.500
92 2300.0 648.031 648.156 0.125
93 2325.0 647.938 648.156 0.219
94 2350.0 647.781 648.156 0.375
95 2375.0 647.688 648.156 0.469
96 2400.0 648.031 648.156 0.125
97 2425.0 647.969 648.156 0.188
98 2450.0 647.875 648.156 0.281
99 2475.0 648.125 648.156 0.031
100 2500.0 648.062 648.156 0.094
101 2525.0 647.969 648.156 0.188
102 2550.0 647.875 648.156 0.281
Text Only
Terrain-derived profile written. Minimum crest raised from 647.66 ft to 648.16 ft.

Add Modified Gates

Python
new_gates = [
    {
        "GateName": "Main Gate",
        "Width": 10.0,
        "Height": 15.0,
        "InvertElevation": 590.0,
        "GateCoefficient": 0.65,
        "NumOpenings": 2,
        "OpeningStations": [5745.0, 5765.0],
    },
    {
        "GateName": "Auxiliary Gate",
        "Width": 5.0,
        "Height": 8.0,
        "InvertElevation": 595.0,
        "GateCoefficient": 0.60,
        "NumOpenings": 1,
        "OpeningStations": [3500.0],
    },
]

GeomLateral.set_connection_gates(
    geom_file,
    TARGET_CONNECTION,
    new_gates,
    create_backup=False,
)

restored_gates = GeomLateral.get_connection_gates(geom_file, TARGET_CONNECTION)
display(restored_gates[gate_cols])

assert len(restored_gates) == 2
for i, expected in enumerate(new_gates):
    actual = restored_gates.iloc[i]
    assert actual["GateName"] == expected["GateName"], f"Gate {i} name mismatch"
    assert actual["Width"] == expected["Width"], f"Gate {i} width mismatch"
    assert actual["Height"] == expected["Height"], f"Gate {i} height mismatch"
    assert np.isclose(actual["InvertElevation"], expected["InvertElevation"], atol=0.01), (
        f"Gate {i} invert mismatch: expected {expected['InvertElevation']}, got {actual['InvertElevation']}"
    )
    assert np.isclose(actual["GateCoefficient"], expected["GateCoefficient"], atol=0.01), (
        f"Gate {i} coefficient mismatch: expected {expected['GateCoefficient']}, got {actual['GateCoefficient']}"
    )
    assert actual["NumOpenings"] == expected["NumOpenings"], f"Gate {i} num openings mismatch"
    actual_stations = [float(s) for s in actual["OpeningStations"]]
    expected_stations = [float(s) for s in expected["OpeningStations"]]
    assert actual_stations == expected_stations, (
        f"Gate {i} stations mismatch: expected {expected_stations}, got {actual_stations}"
    )
print(f"Wrote {len(restored_gates)} gates to '{TARGET_CONNECTION}' — all fields verified.")
GateName Width Height InvertElevation GateCoefficient NumOpenings OpeningStations
0 Main Gate 10.0 15.0 590.0 0.65 2 [5745.0, 5765.0]
1 Auxiliary Gate 5.0 8.0 595.0 0.60 1 [3500.0]
Text Only
Wrote 2 gates to 'Dam' — all fields verified.

Before/After Profile Comparison

Python
fig, axes = plt.subplots(
    3,
    1,
    figsize=(12, 9.4),
    sharex=False,
    gridspec_kw={"height_ratios": [3.7, 2.5, 1.15]},
)
ax = axes[0]
ax.plot(
    terrain_derived_profile["Station"],
    terrain_derived_profile["Elevation"],
    color="0.30",
    linewidth=1.8,
    linestyle="--",
    label="Terrain-derived crest before floor",
    alpha=0.85,
    zorder=5,
)
ax.plot(
    modified_profile["Station"],
    modified_profile["Elevation"],
    color="#b91c1c",
    linewidth=2.2,
    label="Modified crest after +0.5 ft floor",
    zorder=4,
)
ax.axhline(new_min_elev, color="#b91c1c", linestyle=":", linewidth=1.5, label=f"New minimum {new_min_elev:.2f} ft")
ax.scatter(
    modified_profile.loc[terrain_min_idx, "Station"],
    new_min_elev,
    marker="v",
    s=70,
    color="#b91c1c",
    edgecolor="white",
    linewidth=0.8,
    zorder=7,
    label="Raised terrain minimum",
)
ax.set_title(f"{TARGET_CONNECTION}: Terrain-Derived Crest Before And After Improvement")
ax.set_xlabel("")
ax.set_ylabel("Elevation (ft)")
set_profile_ylim(ax, [terrain_derived_profile, modified_profile])
add_gate_station_markers(ax, restored_gates)
dedupe_legend(ax, loc="lower right")

zoom_pad = 50.0
zoom_min = max(float(raised_points["Station"].min()) - zoom_pad, 0.0)
zoom_max = min(float(raised_points["Station"].max()) + zoom_pad, float(modified_profile["Station"].max()))
zoom = modified_profile[(modified_profile["Station"] >= zoom_min) & (modified_profile["Station"] <= zoom_max)]
terrain_zoom = terrain_derived_profile[
    (terrain_derived_profile["Station"] >= zoom_min)
    & (terrain_derived_profile["Station"] <= zoom_max)
]

ax_zoom = axes[1]
ax_zoom.plot(
    terrain_zoom["Station"],
    terrain_zoom["Elevation"],
    color="0.25",
    linewidth=1.8,
    linestyle="--",
    label="Before",
    zorder=5,
)
ax_zoom.plot(
    zoom["Station"],
    zoom["Elevation"],
    color="#b91c1c",
    linewidth=2.2,
    label="After",
    zorder=4,
)
ax_zoom.fill_between(
    zoom["Station"],
    zoom["TerrainElevation"],
    zoom["Elevation"],
    where=zoom["Delta"] > 0,
    color="#fca5a5",
    alpha=0.6,
    label="Filled low spots",
    zorder=3,
)
ax_zoom.axhline(new_min_elev, color="#b91c1c", linestyle=":", linewidth=1.4)
ax_zoom.annotate(
    f"+{CREST_IMPROVEMENT_FT:.1f} ft minimum",
    xy=(float(modified_profile.loc[terrain_min_idx, "Station"]), new_min_elev),
    xytext=(float(modified_profile.loc[terrain_min_idx, "Station"]) + 120.0, new_min_elev + 1.1),
    arrowprops=dict(arrowstyle="->", color="#b91c1c", linewidth=1.0),
    color="#b91c1c",
    fontsize=8,
    bbox=dict(facecolor="white", edgecolor="#fca5a5", alpha=0.85, pad=2),
)
ax_zoom.set_xlim(zoom_min, zoom_max)
ax_zoom.set_title("Focused View Of +0.5 ft Minimum Crest Improvement", pad=10)
ax_zoom.set_xlabel("")
ax_zoom.set_ylabel("Elevation (ft)")
set_profile_ylim(ax_zoom, [terrain_zoom, zoom])
dedupe_legend(ax_zoom, loc="lower right", fontsize=8)

plot_gate_band(axes[2], restored_gates)
axes[2].set_xlim(
    float(terrain_derived_profile["Station"].min()),
    float(terrain_derived_profile["Station"].max()),
)
fig.subplots_adjust(hspace=0.55)
save_figure(fig, "214_connection_authoring_before_after_profile.png")
plt.show()
Text Only
Saved figure: <workspace>\working\connection_authoring\figures\214_connection_authoring_before_after_profile.png

jpeg

Map View Of Modified Connection

Python
fig, ax = plt.subplots(figsize=(11, 7))
plot_context_map(
    ax,
    terrain_raster,
    mesh_areas,
    mesh_cells,
    bc_lines,
    connection_bounds,
    f"{TARGET_CONNECTION}: Plan View With Crest Raise Highlighted",
    show_cells=True,
)

points = modified_profile[["X", "Y"]].to_numpy()
segments = np.stack([points[:-1], points[1:]], axis=1)
segment_delta = np.maximum(
    modified_profile["Delta"].to_numpy()[:-1],
    modified_profile["Delta"].to_numpy()[1:],
)
norm = Normalize(vmin=0.0, vmax=max(CREST_IMPROVEMENT_FT, float(segment_delta.max())))
line_collection = LineCollection(
    segments,
    cmap="inferno",
    norm=norm,
    linewidth=3.4,
    zorder=12,
)
line_collection.set_array(segment_delta)
ax.add_collection(line_collection)

low_pt = modified_profile.loc[terrain_min_idx]
ax.scatter(
    [low_pt["X"]],
    [low_pt["Y"]],
    marker="v",
    s=90,
    color="#b91c1c",
    edgecolor="white",
    linewidth=0.9,
    zorder=14,
    label="Raised low point",
)
ax.text(
    low_pt["X"],
    low_pt["Y"],
    f"STA {low_pt['Station']:.0f}\n+{CREST_IMPROVEMENT_FT:.1f} ft min",
    fontsize=8,
    ha="left",
    va="bottom",
    bbox=dict(facecolor="white", edgecolor="#b91c1c", alpha=0.85, pad=2),
    zorder=15,
)
cbar = fig.colorbar(line_collection, ax=ax, fraction=0.035, pad=0.02)
cbar.set_label("Crest raise above terrain-derived profile (ft)")
dedupe_legend(ax, loc="upper left")
fig.tight_layout()
save_figure(fig, "214_connection_authoring_modified_connection_map.png")
plt.show()
Text Only
Saved figure: <workspace>\working\connection_authoring\figures\214_connection_authoring_modified_connection_map.png

jpeg

Summary

Python
new_coords = GeomLateral.get_connection_line_coords(geom_file, TARGET_CONNECTION)
final_conns = GeomLateral.get_connections(geom_file)

summary = pd.DataFrame([{
    "Connection": TARGET_CONNECTION,
    "Original Weir Coef": 3.82,
    "Modified Weir Coef": 2.60,
    "Original Weir Width": 100,
    "Modified Weir Width": 120,
    "Terrain Source": str(terrain_raster.relative_to(project_path)),
    "Terrain Samples": len(terrain_derived_profile),
    "Minimum Crest Raise": f"+{CREST_IMPROVEMENT_FT:.1f} ft",
    "Points Floored": int((modified_profile["Delta"] > 0).sum()),
    "Original Gates": len(original_gates),
    "Modified Gates": len(restored_gates),
    "Line Points Preserved": len(new_coords) == len(original_coords),
    "Connections Preserved": len(final_conns) == len(connections),
}])
display(summary)

assert len(final_conns) == len(connections)
assert len(new_coords) == len(original_coords)
assert len(restored_gates) == 2
print("Connection authoring workflow complete.")
Connection Original Weir Coef Modified Weir Coef Original Weir Width Modified Weir Width Terrain Source Terrain Samples Minimum Crest Raise Points Floored Original Gates Modified Gates Line Points Preserved Connections Preserved
0 Dam 3.82 2.6 100 120 Terrain\Terrain50.dtm_20ft.tif 298 +0.5 ft 16 1 2 True True
Text Only
Connection authoring workflow complete.

Validate Write Format With HEC-RAS Compute

The file-level assertions above confirm Python can read back what it wrote, but the ultimate validation is running HEC-RAS on the modified geometry and checking that preprocessing and simulation complete without errors. This catches format issues that round-trip checks cannot: malformed fixed-width encoding, invalid coordinate precision, or structural corruption in the connection block.

The demonstration cells above delete and re-create the Dam connection, which removes breach parameters required by Plan 04 (a dam break simulation). To isolate the profile/gate write validation, we extract a fresh project copy and apply set_connection_profile and set_connection_gates with the modified data — then run HEC-RAS directly on that modified geometry without restoring. This proves the API's write format produces geometry HEC-RAS can preprocess and simulate.

Python
import shutil
import re
from ras_commander import init_ras_project, RasCmdr
from ras_commander.hdf import HdfResultsPlan

logging.getLogger("ras_commander").setLevel(logging.INFO)

# Extract a FRESH project copy for compute validation (independent of demo above)
validate_path = WORK_ROOT / "compute_validation"
if validate_path.exists():
    shutil.rmtree(validate_path)
validate_path.mkdir(parents=True, exist_ok=True)

fresh_project = RasExamples.extract_project(PROJECT_NAME, validate_path, suffix="validate")
fresh_geom = fresh_project / GEOM_FILE_NAME

# --------------------------------------------------------------------------
# Write MODIFIED profile and gates to the fresh geometry (no restore).
# HEC-RAS will run on this modified file to prove the write format is valid.
# --------------------------------------------------------------------------

# Create a validation profile: terrain-derived but floored at 651 ft to stay
# above the mesh cell base elevations (~649.6 ft). This profile differs from
# the original (298 points vs 6, terrain-sampled shape vs flat) while remaining
# physically valid for simulation.
VALIDATION_FLOOR_ELEV = 651.0
validation_profile = modified_profile[["Station", "Elevation"]].copy()
validation_profile["Elevation"] = np.maximum(validation_profile["Elevation"], VALIDATION_FLOOR_ELEV)

# Write the validation profile
GeomLateral.set_connection_profile(
    fresh_geom, TARGET_CONNECTION, validation_profile, create_backup=False
)

# Verify profile round-trip with full assertions
readback_profile = GeomLateral.get_connection_profile(fresh_geom, TARGET_CONNECTION)
assert len(readback_profile) == len(validation_profile), (
    f"Profile length mismatch: wrote {len(validation_profile)}, read {len(readback_profile)}"
)
assert np.isclose(readback_profile["Elevation"].min(), VALIDATION_FLOOR_ELEV, atol=0.001), (
    f"Profile min elevation mismatch: expected {VALIDATION_FLOOR_ELEV}, got {readback_profile['Elevation'].min()}"
)
assert np.allclose(
    readback_profile["Elevation"].values,
    validation_profile["Elevation"].values,
    atol=0.01,
), "Profile elevation values diverged beyond 0.01 ft tolerance"
print(f"Profile write verified: {len(readback_profile)} points, "
      f"min elev {readback_profile['Elevation'].min():.3f} ft, "
      f"max elev {readback_profile['Elevation'].max():.3f} ft")

# Write a modified gate that keeps the original name (Plan 04's unsteady file
# has time series for "Gate #1") but uses completely different dimensions,
# position, and coefficient — proving the gate write format is valid.
validation_gates = [
    {
        "GateName": "Gate #1",
        "Width": 12.0,
        "Height": 18.0,
        "InvertElevation": 585.0,
        "GateCoefficient": 0.70,
        "NumOpenings": 2,
        "OpeningStations": [5700.0, 5800.0],
    },
]
GeomLateral.set_connection_gates(fresh_geom, TARGET_CONNECTION, validation_gates, create_backup=False)

# Verify gate round-trip with FULL field assertions
readback_gates = GeomLateral.get_connection_gates(fresh_geom, TARGET_CONNECTION)
assert len(readback_gates) == len(validation_gates), (
    f"Gate count mismatch: wrote {len(validation_gates)}, read {len(readback_gates)}"
)
for i, expected in enumerate(validation_gates):
    actual = readback_gates.iloc[i]
    assert actual["GateName"] == expected["GateName"], f"Gate {i} name mismatch"
    assert actual["Width"] == expected["Width"], f"Gate {i} width mismatch"
    assert actual["Height"] == expected["Height"], f"Gate {i} height mismatch"
    assert np.isclose(actual["InvertElevation"], expected["InvertElevation"], atol=0.01), (
        f"Gate {i} invert mismatch: expected {expected['InvertElevation']}, got {actual['InvertElevation']}"
    )
    assert np.isclose(actual["GateCoefficient"], expected["GateCoefficient"], atol=0.01), (
        f"Gate {i} coefficient mismatch: expected {expected['GateCoefficient']}, got {actual['GateCoefficient']}"
    )
    assert actual["NumOpenings"] == expected["NumOpenings"], f"Gate {i} num openings mismatch"
    actual_stations = [float(s) for s in actual["OpeningStations"]]
    expected_stations = [float(s) for s in expected["OpeningStations"]]
    assert actual_stations == expected_stations, (
        f"Gate {i} stations mismatch: expected {expected_stations}, got {actual_stations}"
    )
print(f"Gate write verified: {len(readback_gates)} gate(s), all fields match")
print(f"  Modified: Width={validation_gates[0]['Width']}, Height={validation_gates[0]['Height']}, "
      f"Invert={validation_gates[0]['InvertElevation']}, Coef={validation_gates[0]['GateCoefficient']}, "
      f"Openings={validation_gates[0]['NumOpenings']}, Stations={validation_gates[0]['OpeningStations']}")

# --------------------------------------------------------------------------
# Run HEC-RAS on the MODIFIED geometry (not the original).
# This is the definitive format validation ��� HEC-RAS parses the written file.
# --------------------------------------------------------------------------

# Initialize and find plan using g13
ras_val = init_ras_project(fresh_project, "7.0")
plans_using_g13 = ras_val.plan_df[ras_val.plan_df["Geom File"] == "13"]
assert not plans_using_g13.empty, "No plan found using geometry file .g13"
plan_number = plans_using_g13.iloc[0]["plan_number"]

# Update Program Version to match installed HEC-RAS (avoids GUI popup)
plan_file = fresh_project / f"BaldEagleDamBrk.p{plan_number.zfill(2)}"
plan_text = plan_file.read_text()
plan_text = re.sub(r"Program Version=.*", f"Program Version={ras_val.ras_version}", plan_text)
plan_file.write_text(plan_text)

geom_text = fresh_geom.read_text()
geom_text = re.sub(r"Program Version=.*", f"Program Version={ras_val.ras_version}", geom_text)
fresh_geom.write_text(geom_text)
print(f"\nUpdated plan {plan_number} and geometry Program Version to {ras_val.ras_version}")

print(f"Running Plan {plan_number} on MODIFIED geometry (298-pt profile floored at "
      f"{VALIDATION_FLOOR_ELEV} ft + modified gate)...")
RasCmdr.compute_plan(plan_number, ras_object=ras_val, force_geompre=True, num_cores=4)

hdf_path = ras_val.plan_df.loc[
    ras_val.plan_df["plan_number"] == plan_number, "HDF_Results_Path"
].iloc[0]
assert hdf_path and Path(hdf_path).exists(), f"HDF file not created for plan {plan_number}"

messages = HdfResultsPlan.get_compute_messages(Path(hdf_path))
lines = messages.splitlines() if messages else []

# Filter for real errors (exclude iteration convergence table headers and metrics)
errors = [l for l in lines if "ERROR" in l.upper()
          and "VOLUME ACCOUNTING" not in l.upper()
          and "ITERATIONS" not in l.upper()
          and "WSEL ERROR" not in l.upper()]
warnings_list = [l for l in lines if "WARNING" in l.upper()]
unstable = [l for l in lines if "UNSTABLE" in l.upper() or "NEGATIVE DEPTH" in l.upper()]

if errors:
    print(f"\nERRORS ({len(errors)}):")
    for e in errors[:10]:
        print(f"  {e}")
    raise RuntimeError(f"HEC-RAS reported {len(errors)} error(s) with modified geometry")

if unstable:
    print(f"\nInstability warnings ({len(unstable)}) — physics-related, not format errors:")
    for w in unstable[:5]:
        print(f"  {w}")

if warnings_list:
    print(f"\nWarnings ({len(warnings_list)}): first 5 shown below")
    for w in warnings_list[:5]:
        print(f"  {w}")

print(f"\n{'='*70}")
print(f"VALIDATION RESULT")
print(f"{'='*70}")
print(f"  Geometry file:    {fresh_geom.name} (MODIFIED — not original)")
print(f"  Profile:          {len(readback_profile)} points (terrain-derived, floored at {VALIDATION_FLOOR_ELEV} ft)")
print(f"  Gate:             Modified width/height/invert/coef/openings/stations")
print(f"  HEC-RAS version:  {ras_val.ras_version}")
print(f"  Plan {plan_number}:          Preprocessing + simulation COMPLETED")
print(f"  Format errors:    {len(errors)}")
print(f"  HDF results:      {Path(hdf_path).name}")
print(f"{'='*70}")
print(f"  set_connection_profile write format: VALIDATED (HEC-RAS accepted)")
print(f"  set_connection_gates write format:   VALIDATED (HEC-RAS accepted)")
print(f"{'='*70}")
Text Only
Profile write verified: 298 points, min elev 651.000 ft, max elev 692.500 ft
Gate write verified: 1 gate(s), all fields match
  Modified: Width=12.0, Height=18.0, Invert=585.0, Coef=0.7, Openings=2, Stations=[5700.0, 5800.0]



Updated plan 04 and geometry Program Version to 7.0
Running Plan 04 on MODIFIED geometry (298-pt profile floored at 651.0 ft + modified gate)...



======================================================================
VALIDATION RESULT
======================================================================
  Geometry file:    BaldEagleDamBrk.g13 (MODIFIED — not original)
  Profile:          298 points (terrain-derived, floored at 651.0 ft)
  Gate:             Modified width/height/invert/coef/openings/stations
  HEC-RAS version:  7.0
  Plan 04:          Preprocessing + simulation COMPLETED
  Format errors:    0
  HDF results:      BaldEagleDamBrk.p04.hdf
======================================================================
  set_connection_profile write format: VALIDATED (HEC-RAS accepted)
  set_connection_gates write format:   VALIDATED (HEC-RAS accepted)
======================================================================

Verify Structure In HDF Results

The compute validation above proves HEC-RAS accepted the modified geometry. As final proof that the Dam connection actively participated in the simulation, extract its time series data from the HDF results using the structure-level HDF extraction APIs.

Python
from ras_commander.hdf import HdfStruc, HdfResultsBreach

hdf_file = Path(hdf_path)

# List all SA/2D connections with time series in the results
connections_in_results = HdfStruc.list_sa2d_connections(hdf_file)
print(f"SA/2D connections in HDF results: {connections_in_results}")
assert TARGET_CONNECTION in connections_in_results, (
    f"'{TARGET_CONNECTION}' not found in HDF results — expected among {connections_in_results}"
)

# Extract structure-level time series for the Dam connection
structure_ts = HdfResultsBreach.get_structure_variables(hdf_file, TARGET_CONNECTION)
print(f"\nTime series shape: {structure_ts.shape}")
print(f"Columns: {list(structure_ts.columns)}")
print(f"Time range: {structure_ts['datetime'].iloc[0]} → {structure_ts['datetime'].iloc[-1]}")

# Summary statistics proving the connection carried flow
peak_total_flow = structure_ts['total_flow'].abs().max()
peak_weir_flow = structure_ts['weir_flow'].abs().max()
max_hw = structure_ts['hw'].max()
max_tw = structure_ts['tw'].max()

summary = pd.DataFrame([{
    "Connection": TARGET_CONNECTION,
    "Time Steps": len(structure_ts),
    "Peak Total Flow (cfs)": f"{peak_total_flow:,.1f}",
    "Peak Weir Flow (cfs)": f"{peak_weir_flow:,.1f}",
    "Max Headwater (ft)": f"{max_hw:.2f}",
    "Max Tailwater (ft)": f"{max_tw:.2f}",
}])
display(summary)

# Show first and last few time steps
print("\nFirst 5 time steps:")
display(structure_ts.head())
print("\nLast 5 time steps:")
display(structure_ts.tail())

assert peak_total_flow > 0, "No flow through connection — structure did not participate"
print(f"\n{'='*70}")
print(f"STRUCTURE VALIDATION: '{TARGET_CONNECTION}' connection carried flow in simulation")
print(f"  Peak total flow: {peak_total_flow:,.1f} cfs")
print(f"  Peak weir flow:  {peak_weir_flow:,.1f} cfs")
print(f"  Max headwater:   {max_hw:.2f} ft")
print(f"  Max tailwater:   {max_tw:.2f} ft")
print(f"{'='*70}")
Text Only
SA/2D connections in HDF results: ['BaldEagleCr Lower Levee', 'BaldEagleCr Middle Levee', 'BaldEagleCr Upper Levee', 'Dam']

Time series shape: (865, 5)
Columns: ['datetime', 'total_flow', 'weir_flow', 'hw', 'tw']
Time range: 1999-01-01 12:00:00 → 1999-01-04 12:00:00
Connection Time Steps Peak Total Flow (cfs) Peak Weir Flow (cfs) Max Headwater (ft) Max Tailwater (ft)
0 Dam 865 352,906.7 357,950.0 673.01 603.40
Text Only
First 5 time steps:
datetime total_flow weir_flow hw tw
0 1999-01-01 12:00:00 2173.012451 0.0 649.999512 588.433594
1 1999-01-01 12:05:00 2263.506836 0.0 649.996765 588.608032
2 1999-01-01 12:10:00 2354.000732 0.0 649.994202 588.728210
3 1999-01-01 12:15:00 2444.494385 0.0 649.991821 588.841736
4 1999-01-01 12:20:00 2534.989502 0.0 649.989685 588.954590
Text Only
Last 5 time steps:
datetime total_flow weir_flow hw tw
860 1999-01-04 11:40:00 66083.179688 61352.835938 662.003601 596.586853
861 1999-01-04 11:45:00 65969.257812 61239.226562 661.993347 596.578186
862 1999-01-04 11:50:00 65854.046875 61124.335938 661.982971 596.569519
863 1999-01-04 11:55:00 65742.648438 61013.257812 661.972656 596.560791
864 1999-01-04 12:00:00 65628.085938 60899.011719 661.962280 596.552002
Text Only
======================================================================
STRUCTURE VALIDATION: 'Dam' connection carried flow in simulation
  Peak total flow: 352,906.7 cfs
  Peak weir flow:  357,950.0 cfs
  Max headwater:   673.01 ft
  Max tailwater:   603.40 ft
======================================================================
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.