Validation and Calibration Reference¶

pywmp.validation and pywmp.calibration provide tools for comparing model output against gauge and spatial observations, and for automated parameter estimation.

HydroMetrics¶

class HydroMetrics(obs, sim)

Computes standard hydrologic performance metrics given observed and simulated time series.

Parameter	Type	Description
`obs`	array-like	Observed discharge or stage series
`sim`	array-like	Simulated series aligned to `obs`

Metrics returned¶

Metric	Symbol	Formula	Reference
Nash-Sutcliffe Efficiency	NSE	1 − Σ(obs−sim)² / Σ(obs−obs_mean)²	Nash & Sutcliffe, 1970
Kling-Gupta Efficiency	KGE	1 − √[(r−1)² + (α−1)² + (β−1)²]	Gupta et al., 2009
KGE revised	KGE'	Modified α component (CV ratio)	Kling et al., 2012
Root Mean Square Error	RMSE	√(mean[(obs−sim)²])	—
Percent Bias	PBIAS	100 × Σ(sim−obs)/Σobs	Moriasi et al., 2007
Peak flow error	ΔQp	(sim_peak − obs_peak)/obs_peak	—
Peak timing error	Δtp	time(sim_peak) − time(obs_peak)	—

KGE decomposes performance into three components:

r — Pearson correlation coefficient
α — variability ratio (σ_sim / σ_obs, or CV_sim / CV_obs for KGE')
β — bias ratio (μ_sim / μ_obs)

Performance ratings (Moriasi et al., 2007)¶

Rating	NSE	KGE	\|PBIAS\|
Very Good	> 0.75	≥ 0.75	< 10%
Good	0.65–0.75	0.65–0.75	10–15%
Satisfactory	0.50–0.65	0.50–0.65	15–25%
Unsatisfactory	< 0.50	< 0.50	> 25%

Methods¶

HydroMetrics.nse()            -> float
HydroMetrics.kge()            -> float
HydroMetrics.kge_prime()      -> float
HydroMetrics.kge_components() -> KGEComponents
HydroMetrics.rmse()           -> float
HydroMetrics.pbias()          -> float
HydroMetrics.summary()        -> dict
HydroMetrics.rating()         -> str   # "Very Good" | "Good" | "Satisfactory" | "Unsatisfactory"

KGEComponents¶

@dataclass
class KGEComponents:
    kge:   float
    r:     float    # correlation
    alpha: float    # variability ratio
    beta:  float    # bias ratio

ModelVsObserved¶

class ModelVsObserved(station, sim_series, variable='Q')

Pairs simulated output with a USGS NWIS gauge record and computes all HydroMetrics.

Parameter	Description
`station`	`GaugeStation` instance (site number + date range)
`sim_series`	`TimeSeries` from `SimResults`
`variable`	`'Q'` (discharge) or `'H'` (stage)

Example — download, pair, and score

from pywmp.validation import ModelVsObserved, GaugeStation, download_streamflow

station = GaugeStation("02301500", start="2023-08-27", end="2023-08-29")
q_obs   = download_streamflow("02301500", start="2023-08-27", end="2023-08-29")

mv = ModelVsObserved(station, sim_results.outflows["Outlet"])
m  = mv.metrics()
print(m.summary())
# NSE 0.81 · KGE 0.78 · |PBIAS| 6.2 % → Very Good
mv.plot("output/validation_plot.png")

Gauge data is downloaded via USGS Water Services REST API using the same HyRiver infrastructure as pywmp.datasets.

SpatialFloodValidation¶

class SpatialFloodValidation(sim_depth, ref_mask, valid_mask=None, depth_threshold=0.1)

Compares a simulated depth raster against a binary reference flood extent.

Parameter	Type	Description
`sim_depth`	ndarray	Simulated max depth (ft or m)
`ref_mask`	ndarray[bool]	Reference flood extent (True = flooded)
`valid_mask`	ndarray[bool]	Optional AOI mask; only cells within AOI are scored
`depth_threshold`	float	Depth above which a cell is considered flooded (default 0.1 ft)

Contingency table metrics¶

Metric	Definition	Reference
CSI (Critical Success Index)	TP / (TP + FP + FN)	Schaefer, 1990
Hit Rate (HR)	TP / (TP + FN)	Wing et al., 2017
False Alarm Ratio (FAR)	FP / (TP + FP)	Wing et al., 2017
F1 Score	2TP / (2TP + FP + FN)	—

Validation against FEMA NFHL AE zones typically yields CSI ≈ 0.25–0.55 for bathtub-style simulations in flat coastal terrain (Wing et al., 2017). For high-resolution (1–10 m) HLL simulations, CSI > 0.6 indicates good agreement with SAR-derived extents (Bates et al., 2010).

Methods¶

SpatialFloodValidation.metrics()              -> dict   # CSI, HR, FAR, F1, TP, FP, FN, TN
SpatialFloodValidation.plot_contingency(ax=None)
SpatialFloodValidation.export_geotiff(path)            # TP/FP/FN/TN raster

Example — compare simulation against FEMA flood zones

from pywmp.validation import SpatialFloodValidation
import geopandas as gpd

fema = gpd.read_file("data/fema_sfha.gpkg")

val = SpatialFloodValidation(
    simulated_depth_tif="output/flood_depth_100yr.tif",
    reference=fema,
    depth_threshold_ft=0.5,
)
m = val.metrics()
print(f"CSI {m['CSI']:.3f}  HR {m['HitRate']:.3f}  FAR {m['FAR']:.3f}")
val.plot_contingency()  # Green = TP, Blue = FN, Red = FP

Suspended sediment rating curve¶

download_ssc(site_no, start, end)            -> TimeSeries
fit_rating_curve(q, ssc)                     -> RatingCurve
estimate_annual_load_rating(q_annual, curve) -> float  # tons/yr

Example — fit SSC rating curve from USGS data

from pywmp.validation import download_ssc, download_streamflow, \
                             fit_rating_curve, estimate_annual_load_rating

ssc = download_ssc("02301500", start="2020-01-01", end="2023-12-31")
q   = download_streamflow("02301500", start="2020-01-01", end="2023-12-31")

curve = fit_rating_curve(q=q, ssc=ssc)
print(f"SSC = {curve.a:.3f} × Q^{curve.b:.3f}")

load = estimate_annual_load_rating(q_annual=q, rating_curve=curve)
print(f"Mean annual load: {load:.1f} tons/yr")

CalibrationEngine¶

class CalibrationEngine(model_fn, observed, params, objective='NSE',
                         method='differential_evolution', **solver_kwargs)

Optimises params by minimising/maximising the objective function computed from model_fn output vs. observed.

Parameter	Description
`model_fn`	Callable that maps `ParameterSet → TimeSeries`; runs the simulation
`observed`	Observed `TimeSeries`
`params`	List of `ParameterBound` instances
`objective`	`'NSE'`, `'KGE'`, `'KGE_prime'`, `'RMSE'`, or `'PBIAS'`
`method`	`'differential_evolution'` (default), `'sce_ua'`, or `'nelder_mead'`

Solver choices¶

Method	Reference	Notes
`differential_evolution`	Storn & Price, 1997	Global; recommended for >3 parameters
`sce_ua`	Duan et al., 1992	Classic hydrology calibration; 3–15 parameters
`nelder_mead`	Nelder & Mead, 1965	Local; use for final refinement near optimum

Methods¶

CalibrationEngine.run(max_iter=500, popsize=15, tol=1e-4) -> CalibrationResults
CalibrationEngine.run_parallel(n_workers=4, ...)           -> CalibrationResults

ParameterBound¶

@dataclass
class ParameterBound(name, lo, hi, initial=None, log_scale=False)

Field	Description
`name`	Parameter name string (must match `model_fn` kwarg)
`lo`, `hi`	Lower and upper bounds
`initial`	Starting point; if `None`, midpoint of [lo, hi]
`log_scale`	If `True`, search on log₁₀ scale (useful for Ks and similar order-of-magnitude parameters)

ParameterSet¶

@dataclass
class ParameterSet:
    values: dict[str, float]   # {name: optimised_value}
    objective_value: float
    n_evaluations: int
    converged: bool

CalibrationResults¶

@dataclass
class CalibrationResults:
    best: ParameterSet
    all_trials: list[ParameterSet]
    convergence_curve: ndarray   # objective per iteration
    runtime_s: float

GaugeStation¶

@dataclass
class GaugeStation(site_no, start, end, variable='Q')

Thin wrapper for a USGS NWIS gauge record. Lazy-loads data on first access via the USGS Water Services REST API.

Download helpers¶

download_streamflow(site_no, start, end) -> TimeSeries
download_stage(site_no, start, end)      -> TimeSeries
download_ssc(site_no, start, end)        -> TimeSeries   # suspended sediment

Full calibration example

from pywmp.calibration import CalibrationEngine, ParameterBound
from pywmp.validation import download_streamflow, HydroMetrics

params = [
    ParameterBound("cn",     lo=70,  hi=98,  initial=85),
    ParameterBound("lag_hr", lo=1.0, hi=8.0, initial=4.0),
    ParameterBound("K_hr",   lo=0.3, hi=3.0, initial=1.0, log_scale=True),
]

def run_model(cn, lag_hr, K_hr):
    # build and run simulation, return TimeSeries
    ...

engine = CalibrationEngine(
    model_fn=run_model,
    observed=download_streamflow("02290000", "2019-01-01", "2020-12-31"),
    params=params,
    objective="KGE",
    method="differential_evolution",
)
results = engine.run(max_iter=400)
print(f"Best KGE: {results.best.objective_value:.3f}")
print(f"Params:   {results.best.values}")

References¶

Bates, P. D., Horritt, M. S., & Fewtrell, T. J. (2010). A simple inertial formulation of the shallow water equations for efficient two-dimensional flood inundation modelling. Journal of Hydrology, 387(1–2), 33–45. DOI ↗

Duan, Q., Sorooshian, S., & Gupta, V. K. (1992). Effective and efficient global optimization for conceptual rainfall-runoff models. Water Resources Research, 28(4), 1015–1031. DOI ↗

Federal Emergency Management Agency. (2023). National Flood Hazard Layer (NFHL). FEMA MSC ↗

Gupta, H. V., Kling, H., Yilmaz, K. K., & Martinez, G. F. (2009). Decomposition of the mean squared error and NSE performance criteria. Journal of Hydrology, 377(1–2), 80–91. DOI ↗

Kling, H., Fuchs, M., & Paulin, M. (2012). Runoff conditions in the upper Danube basin under an ensemble of climate change scenarios. Journal of Hydrology, 424–425, 264–277. DOI ↗

Moriasi, D. N., et al. (2007). Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Transactions of the ASABE, 50(3), 885–900. DOI ↗

Nash, J. E., & Sutcliffe, J. V. (1970). River flow forecasting through conceptual models. Journal of Hydrology, 10(3), 282–290. DOI ↗

Nelder, J. A., & Mead, R. (1965). A simplex method for function minimization. The Computer Journal, 7(4), 308–313. DOI ↗

Schaefer, J. T. (1990). The critical success index as an indicator of warning skill. Weather and Forecasting, 5(4), 570–575. DOI ↗

Storn, R., & Price, K. (1997). Differential evolution — A simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization, 11(4), 341–359. DOI ↗

U.S. Geological Survey. (2023). National Water Information System (NWIS). USGS NWIS ↗

Wing, O. E. J., et al. (2017). Validation of a 30 m resolution flood hazard model of the conterminous United States. Water Resources Research, 53(9), 7968–7986. DOI ↗