Functions and Debugging

Python

CFS_TO_M3S = 0.028316846592

def cfs_to_m3s(discharge_cfs):
    """Convert discharge from cubic feet per second to cubic metres per second."""
    return discharge_cfs * CFS_TO_M3S

q = cfs_to_m3s(1240.0)
print(f"1240 cfs = {q:.2f} m3/s")

Python

def mean_flow_printed(df):
    print(df["discharge_cfs"].mean())    # prints, returns None

def mean_flow(df):
    return df["discharge_cfs"].mean()    # returns a usable value

print("printed version handed back:", mean_flow_printed(stream))
flow = mean_flow(stream)
print(f"returned version: {flow:.1f} cfs = {cfs_to_m3s(flow):.1f} m3/s")

Python

def sample_id(site, lat, lon):
    return f"{site}_{lat:.2f}_{lon:.2f}"

# positional: order is everything
print(sample_id("IOWARIVER", 41.66, -91.53))

# keyword: order stops mattering
print(sample_id(lon=-91.53, site="IOWARIVER", lat=41.66))

# swapped, and Python cannot tell
print(sample_id("IOWARIVER", -91.53, 41.66))

Python

def mean_flow_by(df, date_col="date", flow_col="discharge_cfs", freq="Y"):
    """Mean discharge grouped by period. Returns a Series indexed by period."""
    out = df.copy()
    key = out[date_col].dt.to_period(freq).astype(str)
    return out.groupby(key)[flow_col].mean()

print(mean_flow_by(stream).head(3))          # yearly, the default
print(mean_flow_by(stream, freq="M").head(3))  # monthly, one keyword changed

Python

def add_site_buggy(name, sites=[]):        # one list, shared forever
    sites.append(name)
    return sites

def add_site(name, sites=None):            # a fresh list every call
    sites = [] if sites is None else sites
    sites.append(name)
    return sites

print("buggy:", add_site_buggy("Clear Creek"), add_site_buggy("Rapid Creek"))
print("fixed:", add_site("Clear Creek"), add_site("Rapid Creek"))

Python

FREEZING_C = 0.0                       # global

def days_at_or_below_freezing(temps):
    count = 0                          # local to this call
    for t in temps:
        if t <= FREEZING_C:            # reads the global, cannot rebind it
            count += 1
    return count

n = days_at_or_below_freezing(weather["temp_min_c"])
print(f"{n} of {len(weather)} Iowa City days at or below freezing")

Python

def mean_of(*values):                  # values arrives as a tuple
    return sum(values) / len(values)

print(mean_of(2.1, 3.4, 2.8))
print(mean_of(*stream["discharge_cfs"].head(5)))    # unpack at the call

def plot_flow(df, **style):            # style arrives as a dict
    fig, ax = plt.subplots(figsize=(7, 2.2))
    ax.plot(df["date"], df["discharge_cfs"], **style)   # and unpacks again
    ax.set_ylabel("Discharge (cfs)")

plot_flow(stream.head(365), color="steelblue", linewidth=0.8)

Python

def seven_day_low(df, flow_col="discharge_cfs"):
    """Lowest 7-day rolling mean discharge, in cubic feet per second.

    Parameters
    ----------
    df : daily record, one row per day
    flow_col : name of the discharge column (cfs)
    """
    return df[flow_col].rolling(7).mean().min()

print(f"Lowest 7-day mean: {seven_day_low(stream):.1f} cfs")
help(seven_day_low)

Python

from pathlib import Path

def peak_above(df: "pd.DataFrame", threshold_cfs: float = 10000.0) -> int:
    """Count days whose discharge exceeds a threshold."""
    return int((df["discharge_cfs"] > threshold_cfs).sum())

print("Days above 10000 cfs:", peak_above(stream))
print("Annotations:", peak_above.__annotations__)

Python

def annual_mean_flow_simple(df):
    """Annual mean discharge (cfs) from daily data."""
    data = df.copy()                   # never mutate the caller's table
    data["year"] = data["date"].dt.year
    return data.groupby("year", as_index=False)["discharge_cfs"].mean()

tiny = pd.DataFrame({"date": pd.to_datetime(["2020-01-01", "2020-01-02", "2021-01-01"]),
                     "discharge_cfs": [10.0, 20.0, 40.0]})
got = annual_mean_flow_simple(tiny)
assert got.loc[got["year"] == 2020, "discharge_cfs"].iloc[0] == 15.0   # (10+20)/2
print(got)

Python

from earthcourse.hydro import annual_mean_flow

mine = annual_mean_flow_simple(stream)
package = annual_mean_flow(stream)
check = mine.merge(package, on="year")

assert np.allclose(check["discharge_cfs"], check["mean_discharge_cfs"])
print("My function matches the package function for all",
      len(check), "years.")

Python

Traceback (most recent call last):
  File "load_gauge.py", line 18, in <module>
    mean_flow = summarize("iowa_streamflow_daily_raw.rdb")
  File "load_gauge.py", line 12, in summarize
    values = to_floats(table["42884_00060_00003"])
  File "load_gauge.py", line 7, in to_floats
    return [float(c) for c in cells]
ValueError: could not convert string to float: '14n'

Python

raw = pd.read_csv(RAW_GAUGE, sep="\t", comment="#", dtype=str)

print("row 0      =", raw.iloc[0].to_list())
print("row 1      =", raw.iloc[1].to_list())
print("row 0 type =", type(raw[FLOW_COL].iloc[0]), repr(raw[FLOW_COL].iloc[0]))

try:
    float(raw[FLOW_COL].iloc[0])
except ValueError as error:
    print("Expected ValueError:", error)

Python

coerced = pd.to_numeric(raw[FLOW_COL], errors="coerce")
print("rows            :", len(coerced))
print("NaN after coerce:", int(coerced.isna().sum()))
print(f"mean with coerce: {coerced.mean():.2f} cfs")
print(f"true mean       : {stream['discharge_cfs'].mean():.2f} cfs")

# A naive filter that also fails: '14n' starts with a digit
kept = raw[FLOW_COL].str.match(r"^-?\d", na=False)
print("rows kept by a starts-with-a-digit filter:", int(kept.sum()), "of", len(raw))

Python

def cfs_to_m3s_buggy(discharge_cfs):
    """Convert cfs to m3/s. This version has a bug on purpose."""
    return discharge_cfs * 0.28316846592

print(f"Buggy   : 1000 cfs -> {cfs_to_m3s_buggy(1000.0):.1f} m3/s")
print(f"Correct : 1000 cfs -> {cfs_to_m3s(1000.0):.1f} m3/s")
print()
print("What the data itself says:")
print(stream[["discharge_cfs", "discharge_m3s"]].head(3).to_string(index=False))

Python

import logging
logging.basicConfig(level=logging.DEBUG, format="%(levelname)s %(message)s")

def read_gauge(path, flow_col):
    table = pd.read_csv(path, sep="\t", comment="#", dtype=str)
    logging.debug("rows read = %d", len(table))
    if any(str(v).endswith(("s", "d")) for v in table.iloc[0]):
        logging.warning("dropping USGS type row: %s", table.iloc[0].to_list())
        table = table.iloc[1:].reset_index(drop=True)
    logging.info("returning %d clean rows", len(table))
    return table.assign(discharge_cfs=table[flow_col].astype(float))

Python

def to_floats(cells):
    breakpoint()
    return [float(c) for c in cells]

to_floats(raw["42884_00060_00003"])

Python

(Pdb) p len(cells)
4019
(Pdb) p cells.iloc[0]
'14n'
(Pdb) p cells.iloc[1]
'1050'
(Pdb) p float(cells.iloc[1])
1050.0
(Pdb) p cells.iloc[0].isdigit()
False
(Pdb) q

Python

# tests/test_gauge.py, run with: pytest -q
def test_type_row_is_dropped():
    table = read_gauge(RAW_GAUGE, FLOW_COL)
    assert table["discharge_cfs"].iloc[0] == 1050.0

def test_row_count_matches_processed():
    assert len(read_gauge(RAW_GAUGE, FLOW_COL)) == len(stream)

test_type_row_is_dropped()
test_row_count_matches_processed()
print("2 passed")

Python

from earthcourse.hydro import flow_duration_curve

fdc = flow_duration_curve(stream)
for pct in (10, 50, 90):
    q = np.interp(pct / 100, fdc["exceedance_probability"], fdc["discharge_cfs"])
    print(f"Q{pct:<3} exceeded {pct:>2}% of days: {q:8.1f} cfs")

fig, ax = plt.subplots(figsize=(7, 2.6))
ax.semilogy(fdc["exceedance_probability"] * 100, fdc["discharge_cfs"])
ax.set_xlabel("Percent of days flow is exceeded")
ax.set_ylabel("Discharge (cfs)")
ax.grid(True, alpha=0.3)