Skip to content

Manning's n Calibration with RasCalibrate

RasCalibrate provides parameter estimation workflows for HEC-RAS models, composing RasPermutation, HDF extraction helpers, and statistical metrics into a streamlined calibration API.

This notebook demonstrates two approaches:

  • Example 1 — Grid Search: evaluate every combination of a Manning's n parameter grid, score against observed WSE, and identify the best fit.
  • Example 2 — Nelder-Mead Optimization: gradient-free optimization starting from a perturbed initial guess, converging toward the true parameter values.

What You'll Learn

  • Define CalibrationPoint objects with observed WSE data
  • Use make_mannings_apply_fn() to create a geometry parameter modifier
  • Run RasCalibrate.grid_search() for exhaustive parameter sweeps
  • Run RasCalibrate.optimize() for gradient-free optimization
  • Visualize the calibration objective landscape and convergence history

Test Model

BaldEagleCrkMulti2D — plan 03 (Single 2D Area), ~19,600 mesh cells, 3-day unsteady simulation (~2.2 minutes per run on 4 cores).

Estimated Runtime

  • Setup + baseline run: ~3 minutes
  • Example 1 (grid search, 9 runs, 4 workers): ~7 minutes
  • Example 2 (optimization, 8 iterations): ~18 minutes
  • Total: ~28 minutes
Python
from ras_commander import (
    init_ras_project, RasCmdr, RasExamples,
    RasCalibrate, CalibrationPoint, make_mannings_apply_fn,
)
from ras_commander.hdf import HdfResultsQuery
from pathlib import Path

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import display

plt.style.use("default")
pd.set_option("display.max_columns", 30)
pd.set_option("display.width", 120)
pd.set_option("display.float_format", "{:.4f}".format)

Step 1 — Extract Project and Run Baseline

The baseline run represents the "true" model state. Its output WSE values serve as the synthetic observed data that the calibration tries to reproduce.

In a real calibration, replace these extracted WSE values with gauge-measured stage records.

Python
# Extract a fresh copy of the example project
project_path = RasExamples.extract_project("BaldEagleCrkMulti2D", suffix="calibration")
ras = init_ras_project(project_path, "6.6")

print("Project initialized:", project_path.name)
print("Plans available:")
display(ras.plan_df[["plan_number", "Plan Title", "flow_type"]].reset_index(drop=True))
Python
# Run plan 03 — Single 2D Area (g09, 3-day flood, ~2.2 min on 4 cores)
# This is the baseline: original Manning's n values from the geometry file
print("Running baseline (plan 03, original Manning's n)...")
result = RasCmdr.compute_plan("03", ras_object=ras, force_rerun=True, num_cores=4)
print(f"Baseline: {result}")

# plan_df is automatically refreshed after compute_plan()
baseline_hdf = Path(
    ras.plan_df.loc[ras.plan_df["plan_number"] == "03", "HDF_Results_Path"].values[0]
)
print(f"Results HDF: {baseline_hdf.name} ({baseline_hdf.stat().st_size / 1e6:.1f} MB)")

Step 2 — Extract Observed WSE at Calibration Locations

Three locations along the Bald Eagle Creek 2D mesh simulate gauge stations. HdfResultsQuery.query_points() snaps each (x, y) coordinate to the nearest mesh cell and returns the WSE at the requested time step.

Python
# Representative (x, y) locations along the simulated reach
# These coordinates are in the model's projected coordinate system (ft)
gauge_locations = {
    "Upper Gauge":  (2048000.0, 355500.0),
    "Middle Gauge": (2024000.0, 335250.0),
    "Lower Gauge":  (1988750.0, 307750.0),
}

# Extract WSE at the final time step (peak of the 3-day flood)
observed_raw = HdfResultsQuery.query_points(
    baseline_hdf,
    list(gauge_locations.values()),
    variable="wse",
    time_index=-1,
    ras_object=ras,
)

observed_raw.insert(0, "gauge", list(gauge_locations.keys()))
print("Observed WSE from baseline run (serve as calibration targets):")
display(observed_raw[["gauge", "x", "y", "value", "cell_id", "mesh_name"]])

Step 3 — Define CalibrationPoint Objects

CalibrationPoint pairs an observation with metadata that RasCalibrate uses to extract modeled values from each plan HDF.

  • extraction_method="2d_cell" — locate the nearest 2D mesh cell by (x, y)
  • time_index=-1 — compare at the final (peak) time step
  • metric="rmse" — per-point RMSE contribution to the composite objective
Python
calibration_points = []
for i, (name, (x, y)) in enumerate(gauge_locations.items()):
    observed_wse = float(observed_raw.iloc[i]["value"])
    pt = CalibrationPoint(
        name=name,
        variable="wse",
        extraction_method="2d_cell",
        observed=observed_wse,
        x=x,
        y=y,
        time_index=-1,
        metric="rmse",
        weight=1.0,
    )
    calibration_points.append(pt)
    print(f"{name:>14s} — observed WSE = {observed_wse:.3f} ft")

Step 4 — Define the Manning's n Apply Function

make_mannings_apply_fn() returns a callable that receives a plan path and a pd.Series of parameter values, reads the 2D land cover Manning's n table from the geometry file, updates the named classes, and writes the file back.

The geometry for plan 03 (BaldEagleDamBrk.g09) has 16 NLCD-derived land cover classes. We calibrate two:

Parameter Land Cover Class Baseline n
n_channel Open Water 0.035
n_forest Deciduous Forest 0.100
Python
# Factory returns an apply_fn that modifies the geometry file in-place
apply_fn = make_mannings_apply_fn(
    zone_column_map={
        "n_channel": "Open Water",       # river channel roughness
        "n_forest":  "Deciduous Forest",  # forested floodplain roughness
    },
    path="plaintext",  # reads/writes the plain-text .g## geometry file
)

print("apply_fn ready.")
print("  n_channel → 'Open Water'       (baseline n = 0.035)")
print("  n_forest  → 'Deciduous Forest' (baseline n = 0.100)")

A 3 × 3 Cartesian grid spanning ±25 % of each baseline value.

Text Only
n_channel  : [0.025, 0.035, 0.045]   (baseline = 0.035)
n_forest   : [0.080, 0.100, 0.120]   (baseline = 0.100)
combinations: 9 plans, 4 parallel workers

RasCalibrate.grid_search() internally: 1. Calls RasPermutation.generate_plans() to copy the project and clone plan 03 into 9 variants, each with apply_fn applied for one parameter combination. 2. Executes all variants in parallel with RasPermutation.execute_and_summarize(). 3. Extracts modeled WSE from each HDF and scores against the calibration points.

Expected outcome: the combination (0.035, 0.100) — the original values — produces the lowest RMSE because its simulation exactly reproduces the baseline.

Python
parameters = {
    "n_channel": [0.025, 0.035, 0.045],
    "n_forest":  [0.080, 0.100, 0.120],
}

n_combinations = len(parameters["n_channel"]) * len(parameters["n_forest"])
print(f"Grid search: {n_combinations} combinations, 4 workers")
print(f"Estimated runtime: ~{int(n_combinations / 4 + 0.5) * 3} minutes\n")

results_df = RasCalibrate.grid_search(
    template_plan="03",
    parameters=parameters,
    apply_fn=apply_fn,
    calibration_points=calibration_points,
    metric="rmse",
    suffix="cal",
    max_workers=4,
    num_cores=4,
    ras_object=ras,
)

print("Grid search complete!")
Python
# Display results (best fit at top)
param_cols = ["n_channel", "n_forest"]
score_cols = ["overall_objective", "n_points_scored"]
point_cols = [c for c in results_df.columns if c.startswith("point_") and c.endswith("_objective")]

print("Grid search results (sorted by RMSE — lower is better):")
display(results_df[param_cols + score_cols + point_cols].round(4))

best = results_df.iloc[0]
print(f"\nBest fit: n_channel={best['n_channel']:.3f}, "
      f"n_forest={best['n_forest']:.3f}, "
      f"RMSE={best['overall_objective']:.4f} ft")
Python
# Heatmap of the objective function landscape
pivot = results_df.pivot(
    index="n_forest", columns="n_channel", values="overall_objective"
)

fig, ax = plt.subplots(figsize=(7, 5))
objective_values = pivot.to_numpy(dtype=float)
finite_objectives = objective_values[np.isfinite(objective_values)]
if finite_objectives.size == 0:
    raise ValueError("Grid search produced no finite RMSE values to plot")
max_objective = max(float(finite_objectives.max()), 0.0)
colorbar_vmax = max(max_objective, 1e-6)
im = ax.imshow(objective_values, cmap="YlOrRd_r", aspect="auto",
               vmin=0, vmax=colorbar_vmax)
cbar = plt.colorbar(im, ax=ax)
cbar.set_label("RMSE (ft)")
if max_objective == 0.0:
    cbar.set_ticks([0.0])
    cbar.set_ticklabels(["0.000"])

ax.set_xticks(range(len(pivot.columns)))
ax.set_yticks(range(len(pivot.index)))
ax.set_xticklabels([f"{v:.3f}" for v in pivot.columns])
ax.set_yticklabels([f"{v:.3f}" for v in pivot.index])
ax.set_xlabel("n_channel  (Open Water Manning's n)")
ax.set_ylabel("n_forest  (Deciduous Forest Manning's n)")
ax.set_title("Grid Search Objective Landscape")

# Annotate each cell with its RMSE
for i in range(len(pivot.index)):
    for j in range(len(pivot.columns)):
        val = objective_values[i, j]
        color = "white" if max_objective > 0 and val < max_objective * 0.4 else "black"
        ax.text(j, i, f"{val:.3f}", ha="center", va="center",
                fontsize=9, color=color)

# Highlight best combination
best_col_idx = list(pivot.columns).index(float(best["n_channel"]))
best_row_idx = list(pivot.index).index(float(best["n_forest"]))
ax.add_patch(plt.Rectangle(
    (best_col_idx - 0.5, best_row_idx - 0.5), 1, 1,
    fill=False, edgecolor="lime", linewidth=2.5, label="Best fit"
))
ax.legend(loc="upper right")

plt.tight_layout()
plt.savefig("calibration_grid_search.png", dpi=150, bbox_inches="tight")
plt.show()
print("Figure saved: calibration_grid_search.png")

Example 2 — Nelder-Mead Optimization

The optimizer starts from the midpoint of the parameter bounds (not the true values) and evaluates RasCalibrate.evaluate_single() iteratively.

Each iteration modifies the geometry, runs plan 03, extracts modeled WSE, and returns the RMSE. Nelder-Mead adjusts the simplex toward the minimum.

Note: Full convergence of a Nelder-Mead calibration may require 30–100 model evaluations (~1–4 hours for this model). max_iterations=8 here demonstrates the trajectory without waiting for full convergence.

Python
print("Running Nelder-Mead optimization (max 8 iterations)...")
print("  Starting from midpoint: n_channel=0.040, n_forest=0.100")
print("  True parameter values:  n_channel=0.035, n_forest=0.100\n")

opt_result = RasCalibrate.optimize(
    plan_number="03",
    parameter_bounds={
        "n_channel": (0.020, 0.060),
        "n_forest":  (0.060, 0.140),
    },
    apply_fn=apply_fn,
    calibration_points=calibration_points,
    metric="rmse",
    method="nelder-mead",
    max_iterations=8,
    num_cores=4,
    ras_object=ras,
)

print(f"Status: {opt_result['message']}")
print(f"Iterations completed: {opt_result.get('nit', 'N/A')}")
print(f"\nStarting parameters: {opt_result['starting_parameters']}")
print(f"Best parameters:     {opt_result['best_parameters']}")
print(f"Best RMSE:           {opt_result['best_objective']:.4f} ft")
Python
# Convergence and trajectory visualization
history = opt_result["iteration_history"]  # already a DataFrame

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

# Left: RMSE vs function evaluation
ax = axes[0]
evaluation_nums = range(1, len(history) + 1)
ax.plot(evaluation_nums, history["raw_objective"], "b-o", markersize=6,
        markerfacecolor="steelblue", label="Evaluated RMSE")
ax.axhline(0.0, color="green", linestyle="--", alpha=0.7, label="Perfect fit (RMSE = 0)")
rmse_values = pd.to_numeric(history["raw_objective"], errors="coerce").to_numpy(dtype=float)
finite_rmse = rmse_values[np.isfinite(rmse_values)]
rmse_max = max(float(finite_rmse.max()), 0.0) if finite_rmse.size else 0.0
ax.set_ylim(0.0, max(rmse_max * 1.05, 1e-6))
ax.set_xlabel("Function evaluation")
ax.set_ylabel("RMSE (ft)")
ax.set_title("Optimization Convergence")
ax.legend()
ax.grid(True, alpha=0.3)

# Right: parameter space trajectory
ax2 = axes[1]
n = len(history)
cmap = plt.cm.Blues
colors = [cmap(0.3 + 0.7 * i / max(n - 1, 1)) for i in range(n)]

for i in range(n - 1):
    ax2.plot(
        [history.iloc[i]["n_channel"], history.iloc[i + 1]["n_channel"]],
        [history.iloc[i]["n_forest"],  history.iloc[i + 1]["n_forest"]],
        color="lightsteelblue", linewidth=1, zorder=1,
    )

sc = ax2.scatter(
    history["n_channel"], history["n_forest"],
    c=list(range(n)), cmap="Blues", s=60, zorder=3,
    vmin=0, vmax=max(n - 1, 1),
)
plt.colorbar(sc, ax=ax2, label="Function evaluation")

# Mark the true parameter values
ax2.scatter([0.035], [0.100], color="limegreen", marker="*",
            s=250, zorder=5, label="True values (0.035, 0.100)",
            edgecolors="darkgreen", linewidths=1)

# Shade the search bounds
ax2.set_xlim(0.018, 0.062)
ax2.set_ylim(0.055, 0.145)
ax2.axvline(0.035, color="green", linestyle=":", alpha=0.4)
ax2.axhline(0.100, color="green", linestyle=":", alpha=0.4)
ax2.set_xlabel("n_channel  (Open Water Manning's n)")
ax2.set_ylabel("n_forest  (Deciduous Forest Manning's n)")
ax2.set_title("Parameter Space Trajectory")
ax2.legend(loc="upper right")
ax2.grid(True, alpha=0.3)

plt.suptitle(f"Nelder-Mead Optimization ({len(history)} function evaluations)", fontsize=13, y=1.02)
plt.tight_layout()
plt.savefig("calibration_optimization.png", dpi=150, bbox_inches="tight")
plt.show()
print("Figure saved: calibration_optimization.png")

Summary

What This Notebook Demonstrated

Step Code Purpose
1 RasCmdr.compute_plan("03") Baseline run — creates synthetic observed WSE
2 HdfResultsQuery.query_points() Extract modeled WSE at gauge locations
3 CalibrationPoint(...) Define observation with extraction metadata
4 make_mannings_apply_fn() Factory for land-cover Manning's n modifier
5 RasCalibrate.grid_search() Exhaustive 3×3 sweep — finds RMSE = 0 at true values
6 RasCalibrate.optimize() Nelder-Mead — converges toward true values

Supported Calibration Patterns

Parameter apply functions: - make_mannings_apply_fn() — 2D land cover Manning's n (plain-text or HDF sidecar) - make_infiltration_apply_fn() — Curve Number / infiltration parameters - make_composite_apply_fn() — combine multiple apply_fns

Extraction methods (CalibrationPoint.extraction_method): - "2d_cell" — nearest 2D mesh cell by (x, y) coordinates (used here) - "ref_line" — named reference line in the HDF - "ref_point" — named reference point in the HDF - "1d_xs" — 1D cross section by river/reach/station

Scoring metrics: rmse, nse, kge, pbias, mae

Adapting for Your Model

Python
# Real USGS stage observations
import dataretrieval.nwis as nwis
obs = nwis.get_record(sites="01550000", service="iv",
                       start="1999-01-01", end="1999-01-04")
observed_wse = obs["00065"].max()  # maximum stage (ft above datum)

# Replace synthetic observed values with real gauge data
pt = CalibrationPoint(
    name="USGS 01550000 Bald Eagle Creek",
    variable="wse",
    extraction_method="2d_cell",
    observed=observed_wse,
    x=<gauge_easting_ft>,
    y=<gauge_northing_ft>,
    time_index=-1,
    metric="rmse",
)
CLB Engineering Corporation  ·  LLM Forward Engineering
RAS Commander is a free and open-source project maintained by CLB Engineering Corporation. For agencies and firms seeking to modernize H&H workflows with LLM Forward approaches, contact CLB to partner with the engineers who wrote the automation.