annual streamflow hydrograph for Potomac River at Little Falls

Python Code

"""Hydrology — annual streamflow hydrograph for Potomac River at Little Falls."""


import numpy as np
import pandas as pd
import plotly.graph_objects as go
import io

USGS_URL = (
    "https://waterservices.usgs.gov/nwis/dv/"
    "?format=rdb&sites=01646500&parameterCd=00060"
    "&startDT=2023-01-01&endDT=2023-12-31"
)


def parse_rdb(text):
    """Parse USGS RDB format: skip # comments and the type-spec row."""
    lines = text.splitlines()
    data_lines = []
    header = None
    skip_next = False
    for line in lines:
        if line.startswith("#"):
            continue
        if header is None:
            header = line.strip().split("\t")
            skip_next = True
            continue
        if skip_next:
            # This is the type-spec row (e.g. "5s  10d  ...")
            skip_next = False
            continue
        data_lines.append(line)

    if not data_lines or header is None:
        raise ValueError("No data found in RDB response")

    df = pd.read_csv(
        io.StringIO("\n".join(data_lines)),
        sep="\t", header=None, names=header,
        low_memory=False
    )
    return df


def find_discharge_col(df):
    """Find the discharge column (pattern: XXXXXXXX_00060_00003)."""
    for col in df.columns:
        if "00060" in str(col) and "cd" not in str(col).lower():
            return col
    # fallback: numeric column that isn't agency/site/datetime
    for col in df.columns:
        if col not in ("agency_cd", "site_no", "datetime"):
            try:
                pd.to_numeric(df[col].dropna().iloc[:5])
                return col
            except Exception:
                pass
    return None


def generate():
    try:
        import requests
        print("fetching USGS streamflow data …")
        r = requests.get(USGS_URL, verify=False, timeout=60)
        r.raise_for_status()

        df = parse_rdb(r.text)
        print(f"  parsed {len(df)} rows; columns: {list(df.columns)}")

        # Find date and discharge columns
        date_col = next((c for c in df.columns if c.lower() in ("datetime", "date")), None)
        if date_col is None:
            raise ValueError("No datetime column")

        dis_col = find_discharge_col(df)
        if dis_col is None:
            raise ValueError(f"No discharge column found in: {list(df.columns)}")

        print(f"  using discharge column: {dis_col}")

        df[date_col] = pd.to_datetime(df[date_col], errors="coerce")
        df[dis_col]  = pd.to_numeric(df[dis_col], errors="coerce")
        df = df.dropna(subset=[date_col, dis_col]).sort_values(date_col)

        if len(df) < 50:
            raise ValueError(f"Too few rows: {len(df)}")

        dates = df[date_col]
        discharge = df[dis_col]

    except Exception as e:
        print(f"  USGS fetch failed ({e}); using synthetic hydrograph")
        np.random.seed(7)
        dates = pd.date_range("2023-01-01", "2023-12-31", freq="D")
        n = len(dates)
        # Potomac-like: spring freshet + summer low flow + fall storms
        base = 2000
        t = np.arange(n)
        spring = 8000 * np.exp(-((t - 90)**2) / (2 * 40**2))   # spring peak ~day 90
        summer = -1200 * np.sin(np.pi * t / 365)               # summer low
        fall   = 3000 * np.exp(-((t - 280)**2) / (2 * 25**2))  # fall storm
        noise  = np.random.lognormal(0, 0.3, n) * 200
        discharge = np.maximum(base + spring + summer + fall + noise, 100)

    fig = go.Figure()

    # Fill under curve
    fig.add_trace(go.Scatter(
        x=dates, y=discharge,
        mode="lines",
        line=dict(color=TEAL, width=2),
        fill="tozeroy",
        fillcolor="rgba(82,179,208,0.15)",
        name="Daily mean discharge",
        hovertemplate="%{x|%b %d}: %{y:,.0f} cfs<extra></extra>",
    ))

    # Annotate peak
    if hasattr(discharge, "values"):
        peak_idx = discharge.idxmax()
        peak_date = dates.loc[peak_idx] if hasattr(dates, "loc") else dates[discharge.argmax()]
        peak_val  = discharge.max()
    else:
        peak_idx  = np.argmax(discharge)
        peak_date = dates[peak_idx]
        peak_val  = discharge[peak_idx]

    fig.add_annotation(
        x=peak_date, y=float(peak_val),
        text=f"Peak: {float(peak_val):,.0f} cfs",
        showarrow=True, arrowhead=2,
        ax=40, ay=-40, font=dict(size=11),
    )

    fig.update_layout(
        xaxis=dict(title=""),
        yaxis=dict(title="Discharge (cfs)", type="log"),
        margin=dict(t=50, b=50, l=70, r=60),
    )
    return fig


if __name__ == "__main__":
    generate()

Made with Plotly