# /// script
# requires-python = ">=3.10"
# dependencies = ["mcp>=1.2", "duckdb>=1.3", "pandas>=2.0", "pyarrow>=15", "openpyxl>=3.1",
#                 "matplotlib>=3.8", "lasio>=0.31",
#                 "pillow>=10", "pypdf>=4"]
# ///
"""Hydrohistorian data -- an MCP server over the bulk data releases.

Give an AI assistant (Claude Code, Claude Desktop, any MCP client) the
Colorado water record that hydrohistorian.com/data publishes: daily
streamflow, diversions, reservoirs, water rights, the call record, DWR's
stream network, and the documents a RiverWare build rests on -- and, statewide,
groundwater: well inventories around a point (xlsx/csv/geojson with provenance),
hydrographs and trends, permits, geophysical logs (LAS where DWR has it) with
DWR's formation picks, DWR's Denver Basin tool, cross-sections, and a tracer for
scanned logs whose accuracy is measured against DWR's own digital logs.

Safest install -- download once, check it, run your local copy:

    curl -O https://hydrohistorian.com/static/data/mcp/hydrohistorian_mcp.py
    shasum -a 256 hydrohistorian_mcp.py      # compare with hydrohistorian.com/data#mcp
    claude mcp add hydrohistorian -- uv run /full/path/to/hydrohistorian_mcp.py

(`uv run <url>` also works, but runs whatever the URL serves at each start.)

Everything runs on YOUR machine. The first use of a division downloads its
Parquet files (each checked against the release's SHA-256) into
~/.cache/hydrohistorian/; queries run locally in DuckDB over those files,
with DuckDB's own lock set so a query can read the release and nothing
else. Bulk results go to files you name with `export`, never through the
conversation. With no network, the last downloaded release keeps working.

Environment:
    HYDROHISTORIAN_DATA_URL   release root (default https://hydrohistorian.com/static/data/releases)
                              -- a local directory works too
    HYDROHISTORIAN_CACHE      cache dir (default ~/.cache/hydrohistorian)
    HYDROHISTORIAN_FILES_URL  public base of Hydrohistorian's document mirror (default: none;
                              documents then come straight from DWR)
    HYDROHISTORIAN_DWR_DAILY_CAP  on-demand DWR requests per day from this machine (default 300)

Groundwater tools ask DWR directly only for three things, and only when you ask: its Denver
Basin tool, a document not yet in the mirror, and that document's PDF export. Those requests
go from your machine, one every 2 s, capped per day, and never use DWR's metered REST API.
"""
from __future__ import annotations

import hashlib
import json
import os
import re
import shutil
import tempfile
import threading
import urllib.error
import urllib.request
from pathlib import Path

__version__ = "2"

DATA_URL = os.environ.get("HYDROHISTORIAN_DATA_URL",
                          "https://hydrohistorian.com/static/data/releases").rstrip("/")
CACHE = Path(os.environ.get("HYDROHISTORIAN_CACHE", Path.home() / ".cache" / "hydrohistorian"))
MAX_ROWS = 200
QUERY_SECONDS = 60
UA = {"User-Agent": f"hydrohistorian-mcp/{__version__}"}
_NAME_OK = re.compile(r"^[a-z_][a-z0-9_]*$")
NETWORK_NOTE = ("network.stream_mile is DWR's: miles UPSTREAM from the mouth (or the state line); "
                "downstream_node_id walks toward the outlet.")
GUIDE = (
    "Tables and their key columns (hh_describe has every column and meaning): "
    "gages(source, site_id, name, state, water_district, first_date, last_date); "
    "streamflow_daily(source, site_id, date, flow_cfs, qualifier) -- ALWAYS filter source too: "
    "a DWR gage and its USGS twin (e.g. DWR ANIDURCO / USGS 09361500) are the same river; "
    "diversions_daily(wdid, date, record, value, units) -- filter record='Diversion' (or 'Release') "
    "and units ('CFS' daily rate or 'AF' volume); structures(wdid, name, water_source, stream_mile); "
    "reservoirs(reservoir, hydrodata_site); reservoir_daily(reservoir, date, series, value, units, flag) "
    "-- filter series, and exclude flag='implausible'; water_rights; calls(set_at, released_at, ...); "
    "network(node_id, stream, stream_mile, downstream_node_id); streams; area_capacity; "
    "wells; well_levels; documents.")


class UserError(ValueError):
    """A problem the model can fix: said plainly, returned as {'error': ...}."""


# ---------------------------------------------------------------------------
# Getting the files
# ---------------------------------------------------------------------------

def _read_bytes(rel: str, timeout: float = 60) -> bytes:
    if re.match(r"^https?://", DATA_URL):
        req = urllib.request.Request(f"{DATA_URL}/{rel}", headers=UA)
        with urllib.request.urlopen(req, timeout=timeout) as r:
            return r.read()
    return (Path(DATA_URL) / rel).read_bytes()


def releases() -> list[dict]:
    return json.loads(_read_bytes("index.json", timeout=20))["releases"]


def _rel(division) -> str:
    """Release directory: 'div7' for a water division, 'gw' for statewide groundwater."""
    if str(division).strip().lower() == "gw":
        return "gw"
    try:
        return f"div{int(division)}"
    except (TypeError, ValueError):
        raise UserError(f"division is a water division number (1-7) or 'gw', not {division!r}")


def _cached_manifest(division: int | str) -> dict | None:
    p = CACHE / f"{_rel(division)}.manifest.json"
    try:
        return json.loads(p.read_text(encoding="utf-8"))
    except (OSError, ValueError):
        return None


def manifest(division: int | str) -> tuple[dict, bool]:
    """(manifest, is_live). Offline, the last downloaded release is used."""
    try:
        return json.loads(_read_bytes(f"{_rel(division)}/manifest.json", timeout=20)), True
    except Exception as e:  # noqa: BLE001
        m = _cached_manifest(division)
        if m is not None:
            return m, False
        if "404" in str(e):
            try:
                have = ", ".join(str(r["division"]) for r in releases())
            except Exception:  # noqa: BLE001
                have = "unknown (offline)"
            raise UserError(f"Water Division {division} has no release. Divisions with one: {have}.")
        raise UserError(f"Cannot reach {DATA_URL} and nothing is cached yet: {e}")


def _check_manifest(m: dict) -> None:
    """A manifest names files and tables that get written and queried; a
    tampered one must not be able to escape the cache or inject SQL."""
    for t in m.get("tables", {}):
        if not _NAME_OK.match(t):
            raise UserError(f"refusing release: bad table name {t!r}")
    for f in m.get("files", []):
        if f.get("format") == "parquet" and (Path(f["path"]).name != f["path"] or f["path"].startswith(".")):
            raise UserError(f"refusing release: bad file path {f['path']!r}")


def _local(division: int | str) -> tuple[Path, dict, bool]:
    """(directory of verified Parquet, manifest, is_live)."""
    m, live = manifest(division)
    _check_manifest(m)
    tag = hashlib.sha256(m["built_at"].encode()).hexdigest()[:12]
    d = CACHE / f"{_rel(division)}-{tag}"
    want = sorted(f["path"] for f in m["files"] if f["format"] == "parquet")
    if d.is_dir() and sorted(p.name for p in d.glob("*.parquet")) == want:
        return d, m, live
    if not live:
        raise UserError("offline, and the cached release is incomplete")
    CACHE.mkdir(parents=True, exist_ok=True)
    tmp = Path(tempfile.mkdtemp(prefix=f".{_rel(division)}-", dir=CACHE))
    try:
        for f in m["files"]:
            if f["format"] != "parquet":
                continue
            blob = _read_bytes(f"{_rel(division)}/{f['path']}", timeout=120)
            if hashlib.sha256(blob).hexdigest() != f["sha256"]:
                raise UserError(f"checksum mismatch on {f['path']}; the release may be rebuilding -- retry in a few minutes")
            (tmp / f["path"]).write_bytes(blob)
        try:
            tmp.rename(d)
        except OSError:  # another session finished the same release first
            shutil.rmtree(tmp, ignore_errors=True)
    except BaseException:
        shutil.rmtree(tmp, ignore_errors=True)
        raise
    (CACHE / f"{_rel(division)}.manifest.json").write_text(json.dumps(m), encoding="utf-8")
    for old in CACHE.glob(f"{_rel(division)}-*"):  # earlier releases, best effort
        if old != d:
            shutil.rmtree(old, ignore_errors=True)
    return d, m, live


_CONNS: dict[str, tuple] = {}
_LOCK = threading.Lock()


def _db(division: int | str):
    """Views over the release's Parquet. After setup DuckDB may read ONLY the
    cache directory, and its configuration is locked: a query cannot reach
    the rest of your disk or the network."""
    import duckdb
    with _LOCK:
        key = _rel(division)
        if key in _CONNS:
            return _CONNS[key]
        d, m, live = _local(division)
        con = duckdb.connect()
        for table in m["tables"]:
            p = d / f"{table}.parquet"
            if p.exists():
                con.execute(f"CREATE VIEW {table} AS SELECT * FROM read_parquet('{p.as_posix()}')")
        con.execute("SET memory_limit = '2GB'")
        con.execute("SET threads = 4")
        con.execute("SET allowed_directories = ?", [[d.as_posix() + "/"]])
        con.execute("SET enable_external_access = false")
        con.execute("SET lock_configuration = true")
        _CONNS[key] = (con, m, live)
        return _CONNS[key]


def _check_sql(con, sql: str) -> str:
    """One read-only statement, judged by DuckDB's own parser."""
    import duckdb
    try:
        stmts = con.extract_statements(sql)
    except duckdb.Error as e:
        raise UserError(f"SQL does not parse: {e}")
    if len(stmts) != 1:
        raise UserError("one statement at a time")
    kind = str(stmts[0].type).rsplit(".", 1)[-1]
    if kind not in ("SELECT", "EXPLAIN"):
        raise UserError(f"read-only queries only (SELECT/WITH/DESCRIBE); got {kind}")
    return sql


def _run(division: int | str, sql: str):
    con, m, live = _db(division)
    _check_sql(con, sql)
    cur = con.cursor()
    timer = threading.Timer(QUERY_SECONDS, cur.interrupt)
    timer.start()
    try:
        return cur.execute(sql).df(), m, live
    except Exception as e:  # noqa: BLE001
        if "INTERRUPT" in str(e).upper():
            raise UserError(f"query stopped after {QUERY_SECONDS} s -- narrow it (a join on date alone "
                            "multiplies millions of rows)")
        raise UserError(str(e))
    finally:
        timer.cancel()


def _jsonable(df):
    import pandas as pd
    out = df.copy()
    for c in out.columns:
        s = out[c]
        if pd.api.types.is_datetime64_any_dtype(s):
            t = pd.to_datetime(s)
            out[c] = t.dt.strftime("%Y-%m-%d").where(t.dt.normalize() == t, t.astype(str)).where(t.notna(), None)
    return json.loads(out.to_json(orient="values", date_format="iso", default_handler=str))


def _rows(df, limit: int = MAX_ROWS, stale: bool = False) -> dict:
    total = len(df)
    df = df.head(limit)
    out = {"columns": list(df.columns), "rows": _jsonable(df), "row_count": total, "truncated": total > limit}
    if total > limit:
        out["hint"] = "use hh_export to write the full result to a file"
    if stale:
        out["note"] = "offline: answered from the last downloaded release"
    return out


# ---------------------------------------------------------------------------
# Tool implementations (plain functions, so they test without MCP)
# ---------------------------------------------------------------------------

def list_releases() -> dict:
    return {"releases": releases(), "source": f"{DATA_URL}/index.json"}


def describe(division: int | str = 7) -> dict:
    m, live = manifest(division)
    return {"division": m["division"], "built_at": m["built_at"], "live": live,
            "water_districts": m["water_districts"], "sources": m["sources"],
            "tables": {t: {"description": v["description"], "rows": v["rows"],
                           "date_range": v["date_range"], "columns": v["columns"]}
                       for t, v in m["tables"].items()},
            "caveats": m["caveats"], "network_note": NETWORK_NOTE, "guide": GUIDE}


def find(text: str, division: int | str = 7, limit: int = 25) -> dict:
    """Every word must match (so 'Animas at Durango' finds ANIMAS RIVER AT DURANGO)."""
    con, m, live = _db(division)
    limit = max(1, min(int(limit), MAX_ROWS))
    words = [w for w in re.split(r"\s+", text.strip()) if w and w.lower() not in ("at", "near", "the", "of")]
    if not words:
        raise UserError("give a name or id to search for")
    parts = {
        "gages": ("SELECT 'gage' AS kind, source || ' ' || site_id AS id, name, state, water_district, "
                  "first_date, last_date FROM gages", "coalesce(name,'') || ' ' || site_id || ' ' || source"),
        "structures": ("SELECT 'structure' AS kind, wdid AS id, name, 'CO' AS state, water_district, "
                       "NULL AS first_date, NULL AS last_date FROM structures", "coalesce(name,'') || ' ' || wdid"),
        "reservoirs": ("SELECT 'reservoir' AS kind, hydrodata_site AS id, reservoir AS name, NULL AS state, "
                       "NULL AS water_district, first_date, last_date FROM reservoirs",
                       "reservoir || ' ' || riverware_object"),
        "streams": ("SELECT 'stream' AS kind, gnis_id AS id, name, 'CO' AS state, NULL AS water_district, "
                    "NULL AS first_date, NULL AS last_date FROM streams", "coalesce(name,'') || ' ' || gnis_id"),
        "gw_wells": ("SELECT 'monitoring_well' AS kind, CAST(well_id AS VARCHAR) AS id, name, 'CO' AS state, "
                     "water_district, NULL AS first_date, NULL AS last_date FROM gw_wells",
                     "coalesce(name,'') || ' ' || well_id || ' ' || coalesce(receipt,'')"),
        "gw_log_wells": ("SELECT 'log_well' AS kind, CAST(well_id AS VARCHAR) AS id, name, 'CO' AS state, "
                         "water_district, CAST(log_date AS VARCHAR) AS first_date, NULL AS last_date FROM gw_log_wells",
                         "coalesce(name,'') || ' ' || well_id || ' ' || coalesce(receipt,'') || ' ' || coalesce(permit,'')"),
        "gw_permits": ("SELECT 'permit' AS kind, receipt AS id, coalesce(permit,'') || ' ' || coalesce(category,'') "
                       "AS name, 'CO' AS state, water_district, CAST(date_permit_issued AS VARCHAR) AS first_date, "
                       "NULL AS last_date FROM gw_permits", "receipt || ' ' || coalesce(permit,'')"),
    }
    have = {r[0] for r in con.execute("SHOW TABLES").fetchall()}
    sqls, params = [], []
    for t, (select, hay) in parts.items():
        if t in have:
            sqls.append(select + " WHERE " + " AND ".join([f"({hay}) ILIKE ?"] * len(words)))
            params += [f"%{w}%" for w in words]
    df = con.cursor().execute(f"SELECT * FROM ({' UNION ALL '.join(sqls)}) LIMIT {limit}", params).df()
    return _rows(df, limit, stale=not live)


def query(sql: str, division: int | str = 7) -> dict:
    df, m, live = _run(division, sql)
    return _rows(df, stale=not live)


_EXT = {"csv": (".csv",), "parquet": (".parquet",), "riverware": (".txt", ".rw", ".dat")}


def _safe_out(path: str, fmt: str, overwrite: bool, exts: dict | None = None) -> Path:
    """Exports may not land on a dotfile or inside a dot-directory (~/.zshrc,
    ~/.ssh/...), must carry their format's extension, and do not overwrite
    unless asked -- so a steered model cannot turn export into a file write
    anywhere it likes."""
    out = Path(path).expanduser()
    if not out.is_absolute():
        out = Path.cwd() / out
    out = out.resolve()
    if any(part.startswith(".") for part in out.parts[1:]):
        raise UserError("refusing to write a hidden file or into a hidden directory")
    exts = exts or _EXT
    if out.suffix.lower() not in exts[fmt]:
        raise UserError(f"a {fmt} export must end in {' or '.join(exts[fmt])}")
    if out.exists() and not overwrite:
        raise UserError(f"{out} exists; pass overwrite=true to replace it")
    return out


def export(sql: str, path: str, format: str = "csv", division: int | str = 7,
           units: str | None = None, overwrite: bool = False) -> dict:
    """Write a query's full result to a local file (csv, parquet, riverware)."""
    import pandas as pd
    fmt = format.lower()
    if fmt not in _EXT:
        raise UserError("format is csv, parquet or riverware")
    out = _safe_out(path, fmt, overwrite)
    df, m, live = _run(division, sql)
    out.parent.mkdir(parents=True, exist_ok=True)
    if fmt == "csv":
        df.to_csv(out, index=False)
    elif fmt == "parquet":
        df.to_parquet(out, index=False)
    else:
        if not units:
            raise UserError("riverware export needs units, e.g. cfs, acre-ft or ft")
        if not {"date", "value"} <= set(df.columns):
            raise UserError("riverware export needs columns named date and value, e.g. "
                            "SELECT date, flow_cfs AS value FROM streamflow_daily WHERE source='USGS' AND site_id='09361500'")
        v = df.dropna(subset=["value"]).assign(date=lambda x: pd.to_datetime(x["date"]).dt.normalize())
        if v.empty:
            raise UserError("the query returned no values")
        dup = v["date"].duplicated(keep=False)
        if dup.any():
            raise UserError(f"{int(dup.sum())} rows share a date -- the query mixes series. Filter to one "
                            "(source+site_id, record+units, or reservoir+series).")
        s = v.set_index("date")["value"].sort_index()
        full = s.reindex(pd.date_range(s.index[0], s.index[-1], freq="D"))
        lines = [f"# Written by hydrohistorian-mcp {__version__}; missing days are NaN, never zero.",
                 f"start_date: {full.index[0]:%Y-%m-%d} 24:00", f"end_date: {full.index[-1]:%Y-%m-%d} 24:00",
                 "timestep: 1 DAY", f"units: {units}", "scale: 1"]
        lines += ["NaN" if pd.isna(x) else f"{x:.10g}" for x in full.to_numpy()]
        out.write_text("\n".join(lines) + "\n", encoding="utf-8")
    return {"path": str(out), "rows": int(len(df)), "format": fmt,
            **({"note": "offline: from the last downloaded release"} if not live else {})}


def downstream(node_id: str, division: int | str = 7, max_steps: int = 2000) -> dict:
    con, m, live = _db(division)
    cur = con.cursor()
    path, nid, seen = [], node_id, set()
    while nid and nid not in seen and len(path) < max_steps:
        seen.add(nid)
        r = cur.execute("SELECT node_id, node_type, name, stream, stream_mile, mile_source, downstream_node_id "
                        "FROM network WHERE node_id = ?", [nid]).fetchone()
        if r is None:
            break
        path.append(dict(zip(["node_id", "node_type", "name", "stream", "stream_mile", "mile_source"], r[:6])))
        nid = r[6]
    if not path:
        raise UserError(f"{node_id} is not on the network. Node ids are a WDID ('3000500'), 'DWR ANIDURCO' "
                        "or 'USGS 09361500'; use hh_find to get one. Structures without a DWR stream mile "
                        "are not on it.")
    return {"path": path, "steps": len(path), "note": NETWORK_NOTE}


def documents(category: str | None = None, division: int | str = 7) -> dict:
    m, live = manifest(division)
    docs = m.get("documents") or []
    if category:
        docs = [d for d in docs if category.lower() in str(d.get("category", "")).lower()]
    return {"documents": docs, "count": len(docs)}


# ---------------------------------------------------------------------------
# Groundwater (the statewide 'gw' release) -- the consultant's bread and butter
# ---------------------------------------------------------------------------
#
# Bulk tables come from the release like everything else. Three things are
# fetched from DWR on demand, from YOUR machine, only when a tool asks:
# the Denver Basin determination tool, a document not yet in Hydrohistorian's
# mirror, and that document's PDF export. Those calls are throttled (one
# every 2 s), capped per day, and never touch the metered DWR REST API.

GW = "gw"
EARTH_MI = 3958.8
DWR_TOOLS = "https://dwr.state.co.us/Tools"
WEBLINK = "https://dnrweblink.state.co.us/DWR"
FILES_URL = os.environ.get("HYDROHISTORIAN_FILES_URL", "").rstrip("/")  # the public document mirror
DWR_DAILY_CAP = int(os.environ.get("HYDROHISTORIAN_DWR_DAILY_CAP", "300"))
DWR_INTERVAL_S = 2.0
GAP_DAYS = 3 * 365
_GW_EXT = {"csv": (".csv",), "xlsx": (".xlsx",), "geojson": (".geojson", ".json"),
           "png": (".png",), "las": (".las",), "pdf": (".pdf",)}
DB_TOOL_AQUIFERS = ["Upper Dawson", "Lower Dawson", "Denver", "Upper Arapahoe", "Lower Arapahoe",
                    "Laramie-Fox Hills"]
# DWR pick names -> the tool's aquifer. Coarse picks (DAWSON, ARAPAHOE) span
# two of the tool's aquifers and are compared to neither -- shown, not merged.
COARSE_PICK = {"Upper Dawson": "DAWSON", "Lower Dawson": "DAWSON", "Upper Arapahoe": "ARAPAHOE",
               "Lower Arapahoe": "ARAPAHOE"}
PICK_TO_TOOL = {"UPPER DAWSON": "Upper Dawson", "LOWER DAWSON": "Lower Dawson", "DENVER": "Denver",
                "UPPER ARAPAHOE": "Upper Arapahoe", "LOWER ARAPAHOE": "Lower Arapahoe",
                "LARAMIE-FOX HILLS": "Laramie-Fox Hills"}
DENVER_TOOL_NOTE = (
    "Output of DWR's Denver Basin location tool (dwr.state.co.us/Tools/DenverLocation), returned as "
    "DWR computed it. It is DWR's screening estimate from its aquifer surfaces, not a determination "
    "of water availability by the State Engineer and not a substitute for one. 'annual_approp_af' is "
    "DWR's 'Annual Approp. (AF)' column for the acreage given; 'status' is DWR's NT "
    "(nontributary) / NNT (not-nontributary) code.")


def _num_or_none(v):
    try:
        x = float(str(v).strip())
        return None if x != x else x
    except (TypeError, ValueError):
        return None


def _check_point(lat, lon) -> tuple[float, float]:
    lat, lon = float(lat), float(lon)
    if not (36.8 <= lat <= 41.2 and -109.2 <= lon <= -101.9):
        raise UserError(f"({lat}, {lon}) is not in Colorado; give decimal degrees, longitude negative")
    return lat, lon


def _dist_sql(lat: float, lon: float, lat_col="latitude", lon_col="longitude") -> str:
    return (f"{EARTH_MI}*2*asin(sqrt(pow(sin(radians({lat_col}-{lat})/2),2)+cos(radians({lat}))"
            f"*cos(radians({lat_col}))*pow(sin(radians({lon_col}-({lon}))/2),2)))")


def _bbox_sql(lat: float, lon: float, r: float, lat_col="latitude", lon_col="longitude") -> str:
    import math
    dlat = r / 69.0
    dlon = r / max(1e-6, 69.0 * math.cos(math.radians(lat)))
    return (f"{lat_col} BETWEEN {lat - dlat} AND {lat + dlat} AND "
            f"{lon_col} BETWEEN {lon - dlon} AND {lon + dlon}")


def _gw():
    con, m, live = _db(GW)
    return con.cursor(), m, live


def _near(cur, table: str, select: str, lat, lon, r, extra_where: str = "", params=(), order="distance_mi"):
    sql = (f"SELECT round({_dist_sql(lat, lon)}, 3) AS distance_mi, {select} FROM {table} "
           f"WHERE {_bbox_sql(lat, lon, r)} {extra_where}")
    df = cur.execute(f"SELECT * FROM ({sql}) WHERE distance_mi <= {float(r)} ORDER BY {order}",
                     list(params)).df()
    return df


PERMIT_COLS = ("receipt, permit, status, category, uses, aquifers AS aquifer_label_permit, depth_total_ft, "
               "top_perf_ft, bottom_perf_ft, static_water_level_ft, static_water_level_date, "
               "pump_test_yield_gpm, date_well_completed, date_well_plugged, driller, county, "
               "designated_basin, pm, township, range, section, q40, q160, location_accuracy, "
               "latitude, longitude, dwr_url")
MONITOR_COLS = ("well_id, name, receipt, wdid, regime, aquifer1 AS aquifer_label, well_depth_ft, top_perf_ft, "
                "bottom_perf_ft, elevation_ft, designated_basin, latitude, longitude")
LOG_COLS = ("well_id, name, receipt, permit, log_date, log_type_raw, curve_families, log_top_ft, log_bottom_ft, "
            "log_swl_ft, aquifer_picks, well_depth_ft, ecmc_id, ecmc_url, latitude, longitude, dwr_url")


def _inventory(lat, lon, radius_mi, aquifer=None, kinds=("permits", "monitoring", "logs")):
    lat, lon = _check_point(lat, lon)
    r = float(radius_mi)
    if not 0 < r <= 10:
        raise UserError("radius_mi is between 0 and 10 (use hh_query for anything wider)")
    cur, m, live = _gw()
    out = {}
    aq = f"%{aquifer}%" if aquifer else None
    if "permits" in kinds:
        out["permits"] = _near(cur, "gw_permits", PERMIT_COLS, lat, lon, r,
                               "AND aquifers ILIKE ?" if aq else "", [aq] if aq else [])
    if "monitoring" in kinds:
        df = _near(cur, "gw_wells", MONITOR_COLS, lat, lon, r,
                   "AND (aquifer1 ILIKE ? OR aquifer2 ILIKE ? OR regime ILIKE ?)" if aq else "",
                   [aq, aq, aq] if aq else [])
        if not df.empty:
            ids = ",".join(str(int(i)) for i in df["well_id"])
            last = cur.execute(
                f"SELECT well_id, count(*) AS n_levels, min(date) AS first_level, max(date) AS last_level, "
                f"arg_max(depth_to_water_ft, date) AS last_depth_to_water_ft FROM gw_levels "
                f"WHERE well_id IN ({ids}) GROUP BY well_id").df()
            tr = cur.execute(f"SELECT well_id, trend, slope_ft_per_decade, p_value FROM gw_trends "
                             f"WHERE well_id IN ({ids})").df()
            df = df.merge(last, on="well_id", how="left").merge(tr, on="well_id", how="left")
        out["monitoring_wells"] = df
    if "logs" in kinds:
        df = _near(cur, "gw_log_wells", LOG_COLS, lat, lon, r,
                   "AND aquifer_picks ILIKE ?" if aq else "", [aq] if aq else [])
        if not df.empty:
            ids = ",".join(str(int(i)) for i in df["well_id"])
            docs = cur.execute(
                f"SELECT well_id, count(*) AS n_documents, bool_or(doc_type = 'Log LAS') AS has_las "
                f"FROM gw_documents WHERE well_id IN ({ids}) GROUP BY well_id").df()
            df = df.merge(docs, on="well_id", how="left")
        out["log_wells"] = df
    return out, m, live


def wells_near(lat, lon, radius_mi: float = 1.0, aquifer: str | None = None, limit: int = 50) -> dict:
    frames, m, live = _inventory(lat, lon, radius_mi, aquifer)
    limit = max(1, min(int(limit), MAX_ROWS))
    out = {"center": [float(lat), float(lon)], "radius_mi": float(radius_mi), "release_built_at": m["built_at"]}
    for k, df in frames.items():
        out[k] = _rows(df, limit, stale=not live)
    out["notes"] = [
        "permits: one row per receipt (an application, not necessarily a drilled well); group by permit. "
        "Depths, yields and static levels are the driller's reports.",
        "aquifer_label_* are labels as filed; formation depths come from log-well picks (gw_formation_tops).",
        "gw_export_inventory writes all three tables with a provenance sheet.",
    ]
    return out


def _provenance(m: dict, extra: dict | None = None) -> list[tuple[str, str]]:
    rows = [("Source", "Hydrohistorian statewide groundwater release (public DWR records)"),
            ("Release built", m["built_at"]), ("Release manifest", f"{DATA_URL}/{GW}/manifest.json"),
            ("Written by", f"hydrohistorian-mcp {__version__}")]
    rows += [(f"{f['path']} sha256", f["sha256"]) for f in m["files"] if f["format"] == "parquet"]
    rows += [(k, str(v)) for k, v in (extra or {}).items()]
    rows += [("Caveat", c) for c in m.get("caveats", [])]
    return rows


def export_inventory(lat, lon, radius_mi: float, path: str, format: str = "xlsx",
                     aquifer: str | None = None, overwrite: bool = False) -> dict:
    import pandas as pd
    fmt = format.lower()
    if fmt not in ("xlsx", "csv", "geojson"):
        raise UserError("format is xlsx, csv or geojson")
    out = _safe_out(path, fmt, overwrite, exts=_GW_EXT)
    frames, m, live = _inventory(lat, lon, radius_mi, aquifer)
    prov = _provenance(m, {"Center (lat, lon)": f"{lat}, {lon}", "Radius (mi)": radius_mi,
                           "Aquifer filter": aquifer or "none"})
    col_docs = {t: m["tables"][t]["columns"] for t in ("gw_permits", "gw_wells", "gw_log_wells")}
    out.parent.mkdir(parents=True, exist_ok=True)
    written = [str(out)]
    if fmt == "xlsx":
        names = {"permits": "Permits", "monitoring_wells": "Monitoring wells", "log_wells": "Log wells"}
        with pd.ExcelWriter(out, engine="openpyxl") as xw:
            pd.DataFrame(prov, columns=["Item", "Value"]).to_excel(xw, sheet_name="README", index=False)
            for k, df in frames.items():
                df.to_excel(xw, sheet_name=names[k], index=False)
            pd.DataFrame([(t, c, d) for t, cols in col_docs.items() for c, d in cols.items()],
                         columns=["Table", "Column", "Meaning"]).to_excel(xw, sheet_name="Columns", index=False)
    elif fmt == "csv":
        written = []
        for k, df in frames.items():
            p = out.with_name(f"{out.stem}_{k}.csv")
            if p.exists() and not overwrite:
                raise UserError(f"{p} exists; pass overwrite=true")
            df.to_csv(p, index=False)
            written.append(str(p))
        readme = out.with_name(f"{out.stem}_README.txt")
        readme.write_text("\n".join(f"{a}: {b}" for a, b in prov) + "\n", encoding="utf-8")
        written.append(str(readme))
    else:
        feats = []
        for k, df in frames.items():
            for rec in json.loads(df.to_json(orient="records", date_format="iso", default_handler=str)):
                if rec.get("latitude") is None or rec.get("longitude") is None:
                    continue
                feats.append({"type": "Feature", "properties": {"kind": k, **rec},
                              "geometry": {"type": "Point", "coordinates": [rec["longitude"], rec["latitude"]]}})
        out.write_text(json.dumps({"type": "FeatureCollection", "features": feats,
                                   "metadata": dict(prov[:4]) | {"caveats": m.get("caveats", [])}}),
                       encoding="utf-8")
    return {"written": written, "counts": {k: int(len(v)) for k, v in frames.items()},
            **({"note": "offline: from the last downloaded release"} if not live else {})}


def hydrograph(well_ids: list, csv_path: str | None = None, png_path: str | None = None,
               overwrite: bool = False) -> dict:
    ids = sorted({int(i) for i in (well_ids if isinstance(well_ids, (list, tuple)) else [well_ids])})
    if not 1 <= len(ids) <= 20:
        raise UserError("give 1 to 20 well_ids (gw_wells_near lists them)")
    cur, m, live = _gw()
    idl = ",".join(map(str, ids))
    wells = cur.execute(f"SELECT well_id, name, regime, setting, aquifer1, well_depth_ft, elevation_ft "
                        f"FROM gw_wells WHERE well_id IN ({idl})").df()
    if wells.empty:
        raise UserError("none of those ids is a DWR monitoring well in the release")
    levels = cur.execute(f"SELECT well_id, date, depth_to_water_ft, water_level_elevation_ft FROM gw_levels "
                         f"WHERE well_id IN ({idl}) ORDER BY well_id, date").df()
    trends = cur.execute(f"SELECT * FROM gw_trends WHERE well_id IN ({idl})").df()
    stats = levels.groupby("well_id").agg(n=("date", "size"), first=("date", "min"), last=("date", "max"),
                                          min_depth_ft=("depth_to_water_ft", "min"),
                                          max_depth_ft=("depth_to_water_ft", "max")).reset_index()
    summary = wells.merge(stats, on="well_id", how="left").merge(
        trends[["well_id", "trend", "slope_ft_per_decade", "p_value", "span_years"]], on="well_id", how="left")
    res = {"wells": _rows(summary, stale=not live), "written": [],
           "trend_rule": "Mann-Kendall on annual means; a trend is reported only with >=20 readings over "
                         ">=10 years (the site's own gate). Bedrock declines are storage loss, not recharge."}
    if csv_path:
        p = _safe_out(csv_path, "csv", overwrite, exts=_GW_EXT)
        p.parent.mkdir(parents=True, exist_ok=True)
        levels.to_csv(p, index=False)
        res["written"].append(str(p))
    if png_path:
        import matplotlib
        matplotlib.use("Agg")
        import matplotlib.pyplot as plt
        p = _safe_out(png_path, "png", overwrite, exts=_GW_EXT)
        fig, ax = plt.subplots(figsize=(10, 5.5), dpi=150)
        names = dict(zip(wells["well_id"], wells["name"]))
        for wid, g in levels.groupby("well_id"):
            t = pd_dates(g["date"]).reset_index(drop=True)
            y = g["depth_to_water_ft"].reset_index(drop=True).astype(float)
            tl, yl = [], []  # a NaN inside each long gap: no line across years without a reading
            for i in range(len(t)):
                if i and (t[i] - t[i - 1]).days > GAP_DAYS:
                    tl.append(t[i - 1] + (t[i] - t[i - 1]) / 2)
                    yl.append(float("nan"))
                tl.append(t[i])
                yl.append(y[i])
            line, = ax.plot(tl, yl, lw=0.8, label=f"{names.get(wid) or wid} ({wid})")
            ax.plot(t, y, ls="none", marker=".", ms=3, color=line.get_color())
        ax.invert_yaxis()
        ax.text(0.01, 0.98, f"Lines break where readings are more than {GAP_DAYS // 365} years apart",
                transform=ax.transAxes, fontsize=7, va="top", color="0.35")
        ax.set_ylabel("Depth to water, ft below land surface")
        ax.set_title("DWR monitoring wells: depth to water")
        ax.grid(alpha=0.3)
        ax.legend(fontsize=7, loc="best")
        fig.text(0.01, 0.01, f"Source: Colorado DWR via Hydrohistorian release {m['built_at']}", fontsize=6)
        fig.tight_layout()
        p.parent.mkdir(parents=True, exist_ok=True)
        fig.savefig(p)
        plt.close(fig)
        res["written"].append(str(p))
        res["png_alt_text"] = (f"Line chart of depth to water (feet, increasing downward) over time for "
                               f"{len(ids)} DWR monitoring well(s): " + "; ".join(
                                   (f"{r.name or r.well_id}: {int(r.n)} readings {str(r.first)[:10]} to {str(r.last)[:10]}, "
                                    f"{r.min_depth_ft}-{r.max_depth_ft} ft" if r.n == r.n and r.n else
                                    f"{r.name or r.well_id}: no readings")
                                   for r in summary.itertuples()) + ". The CSV holds every value.")
    return res


def pd_dates(s):
    import pandas as pd
    return pd.to_datetime(s)


def permit_fn(receipt: str | None = None, permit: str | None = None) -> dict:
    if not receipt and not permit:
        raise UserError("give receipt or permit (e.g. permit='22189-F')")
    cur, m, live = _gw()
    if receipt:
        df = cur.execute("SELECT * FROM gw_permits WHERE receipt = ?", [str(receipt).strip()]).df()
    else:
        df = cur.execute("SELECT * FROM gw_permits WHERE permit = ? OR permit ILIKE ?",
                         [str(permit).strip(), f"{str(permit).strip()}%"]).df()
    receipts = [r for r in df["receipt"].dropna().unique()] if not df.empty else (
        [str(receipt).strip()] if receipt else [])
    docs = cur.execute(
        "SELECT DISTINCT d.* FROM gw_documents d LEFT JOIN gw_log_wells w ON d.well_id = w.well_id "
        "WHERE d.receipt IN (SELECT unnest(?)) OR w.receipt IN (SELECT unnest(?)) ORDER BY created",
        [receipts, receipts]).df() if receipts else None
    out = {"permits": _rows(df, stale=not live)}
    if docs is not None:
        out["documents"] = _rows(docs)
    if df.empty:
        out["note"] = "not in the release's permit file (it covers what DWR's permit API returns)"
    return out


def logs_near(lat, lon, radius_mi: float = 2.0, family: str | None = None, limit: int = 100) -> dict:
    frames, m, live = _inventory(lat, lon, radius_mi, kinds=("logs",))
    df = frames["log_wells"]
    if family and not df.empty:
        df = df[df["curve_families"].fillna("").str.contains(family.lower())]
    cur, _, _ = _gw()
    if not df.empty:
        ids = ",".join(str(int(i)) for i in df["well_id"])
        picks = cur.execute(f"SELECT well_id, string_agg(aquifer || ' ' || coalesce(cast(top_depth_ft AS INT)::VARCHAR,'?') "
                            f"|| '-' || coalesce(cast(base_depth_ft AS INT)::VARCHAR,'?'), '; ' ORDER BY top_depth_ft) "
                            f"AS picks_ft FROM gw_picks WHERE well_id IN ({ids}) GROUP BY well_id").df()
        df = df.merge(picks, on="well_id", how="left")
        if _has_table(cur, "gw_digitized"):
            dig = cur.execute(f"SELECT well_id, count(*) AS digitized_curves_pass FROM gw_digitized "
                              f"WHERE gate = 'pass' AND well_id IN ({ids}) GROUP BY well_id").df()
            df = df.merge(dig, on="well_id", how="left")
    return {"log_wells": _rows(df, max(1, min(int(limit), MAX_ROWS)), stale=not live),
            "notes": ["log_type_raw is DWR's free-text code (no published legend); curve_families decodes only "
                      "spelled-out tokens, so a family filter misses logs filed as bare letters (R, G, C...).",
                      "has_las = DWR holds a digital LAS file (gw_documents / gw_fetch_document); older logs are scans.",
                      "digitized_curves_pass = curves auto-digitized from the scan that passed the confidence gate "
                      "(gw_digitized_curves / gw_fetch_digitized).",
                      "Oil and gas wells (ecmc_id) keep their logs at ECMC: follow ecmc_url."]}


def _has_table(cur, table: str) -> bool:
    return bool(cur.execute("SELECT count(*) FROM information_schema.tables WHERE table_name = ?", [table]).fetchone()[0])


def digitized_curves(well_id: int | None = None, lat: float | None = None, lon: float | None = None,
                     radius_mi: float = 2.0, gate: str = "pass", limit: int = 200) -> dict:
    """Curves digitized automatically from the mirrored log scans (gw_digitized), for one log well
    or every log well near a point. gate='pass' (default) keeps only curves that cleared the
    confidence gate; 'fail' shows the rest with their reasons; 'all' both."""
    cur, m, live = _gw()
    if not _has_table(cur, "gw_digitized"):
        raise UserError("this release has no gw_digitized table yet (built before the batch digitizer ran)")
    if gate not in ("pass", "fail", "all"):
        raise UserError("gate is pass, fail or all")
    where, params = [], []
    if well_id is not None:
        where.append("d.well_id = ?"); params.append(int(well_id))
    elif lat is not None and lon is not None:
        w = _near(cur, "gw_log_wells", "well_id", lat, lon, radius_mi)
        if w.empty:
            return {"curves": [], "count": 0, "note": f"no log wells within {radius_mi} mi"}
        where.append("d.well_id IN (" + ",".join(str(int(i)) for i in w["well_id"]) + ")")
    else:
        raise UserError("give well_id, or lat and lon")
    if gate != "all":
        where.append("d.gate = ?"); params.append(gate)
    df = cur.execute("SELECT d.doc_id, d.well_id, d.page, d.curve, d.mnemonic, d.unit, d.scale, d.weight, d.v_left, "
                     "d.v_right, d.depth_top_ft, d.depth_base_ft, d.coverage, d.ambiguous, d.n_wraps, d.labels, "
                     "d.max_resid_ft, d.gate, d.reasons, d.las, d.version FROM gw_digitized d WHERE "
                     + " AND ".join(where) + " ORDER BY d.well_id, d.doc_id, d.page, d.curve", params).df()
    return {"curves": _rows(df, max(1, min(int(limit), MAX_ROWS)), stale=not live),
            "count": int(len(df)),
            "notes": ["DIGITIZED FROM RASTER SCANS by hydrohistorian's automatic digitizer -- screening data, not "
                      "the original log. Keep the scan (gw_documents / gw_fetch_document) beside any value you quote.",
                      "gate=pass means: >=4 OCR'd depth labels on one straight axis (worst residual <= 1 ft), "
                      "coverage >= 40%, <= 5% ambiguous rows, scale numbers read at >= 0.80 confidence, and a "
                      "recognized curve name; SP never passes (its printed scale is a span); on Phoenix/Pioneer prints only "
                      "the solid gamma passes. On 402 scan+LAS wells (3,704 pages) 579 curves pass and 7% of those score below r 0.8 "
                      "against DWR's LAS; the benchmark ships "
                      "with the release (digitizer-benchmark.md).",
                      "A printed SP scale is a span, not absolute millivolts: SP values carry an unknown offset.",
                      "gw_fetch_digitized saves the LAS (one file per scan page; every curve of that page is in it)."]}


def fetch_digitized(doc_id: int, page: int, path: str, overwrite: bool = False) -> dict:
    cur, m, live = _gw()
    if not _has_table(cur, "gw_digitized"):
        raise UserError("this release has no gw_digitized table yet")
    df = cur.execute("SELECT DISTINCT las FROM gw_digitized WHERE doc_id = ? AND page = ?", [int(doc_id), int(page)]).df()
    if df.empty:
        raise UserError(f"no digitized LAS for doc {doc_id} page {page} (gw_digitized_curves lists what exists)")
    if not FILES_URL:
        raise UserError("HYDROHISTORIAN_FILES_URL is not set: the digitized LAS files live in the document mirror; "
                        "set it to the mirror's public base URL (hydrohistorian.com/data says it)")
    out = _safe_out(path, "las", overwrite, exts=_GW_EXT)
    req = urllib.request.Request(f"{FILES_URL}/{df['las'].iloc[0]}", headers=UA)
    with urllib.request.urlopen(req, timeout=60) as r:
        data = r.read()
    out.parent.mkdir(parents=True, exist_ok=True)
    out.write_bytes(data)
    return {"written": str(out), "bytes": len(data), "digitized_from_raster": True,
            "note": "Screening data; the LAS's ~Other block carries the calibration and each curve's gate verdict."}


# --- Net sand from a gamma curve ------------------------------------------------------------
#
# The step a modeler takes after the log: a shale index from gamma, a cutoff, beds, and
# net sand per aquifer interval. Conventions, not measurements, so the tool prints every
# one of them, and where DWR has published its own net-sand pick for the same well and
# aquifer (gw_picks.net_sand_ft, a geologist's reading of the same log) it is shown beside
# ours. The cutoff and percentiles below are the ones that matched DWR's picks best over
# 116 wells with both a LAS and a pick (docs/research/net-sand-benchmark.md).

GAMMA_MNEMONICS = ("GR", "NATURAL_GAMMA", "NATURALGAMMA", "GAMMA", "GAMMARAY", "GAMMA_RAY", "NATURAL GAMMA", "SGR", "CGR")
NET_SAND_DEFAULTS = {"cutoff": 0.4, "min_bed_ft": 1.0, "p_clean": 10, "p_shale": 90}


def _gamma_curve(las, name: str | None = None):
    names = [c.mnemonic for c in las.curves]
    if name:
        if name not in names:
            raise UserError(f"{name!r} is not a curve in this LAS ({', '.join(names)})")
        return name
    up = {n.upper().replace("-", "_"): n for n in names}
    for m in GAMMA_MNEMONICS:
        if m in up:
            return up[m]
    for n in names:
        if "GAMMA" in n.upper() or n.upper().startswith("GR"):
            return n
    raise UserError(f"no gamma curve found ({', '.join(names)}); name it with gr_curve")


def net_sand_from_arrays(depth, gr, cutoff: float, min_bed_ft: float, p_clean: float, p_shale: float,
                         intervals: list[tuple[float, float, str]] | None = None) -> dict:
    """Linear gamma index: Vsh = (GR - GR_clean) / (GR_shale - GR_clean), clean and shale lines
    at the p_clean / p_shale percentiles of the gamma over the log (per interval when given);
    sand where Vsh < cutoff; beds thinner than min_bed_ft merged into their neighbours."""
    import numpy as np
    depth = np.asarray(depth, float); gr = np.asarray(gr, float)
    ok = ~np.isnan(gr) & ~np.isnan(depth)
    depth, gr = depth[ok], gr[ok]
    if len(depth) < 10:
        raise UserError("fewer than 10 gamma samples")
    order = np.argsort(depth); depth, gr = depth[order], gr[order]
    step = float(np.median(np.diff(depth))) if len(depth) > 1 else 0.5

    def beds_for(lo, hi):
        k = (depth >= lo) & (depth <= hi)
        if k.sum() < 10:
            return None
        d, g = depth[k], gr[k]
        g_clean, g_shale = np.percentile(g, p_clean), np.percentile(g, p_shale)
        if g_shale - g_clean < 1e-6:
            return None
        vsh = np.clip((g - g_clean) / (g_shale - g_clean), 0, 1)
        sand = vsh < cutoff
        # runs -> beds, then merge thin beds into the neighbour they sit between
        beds = []
        start = 0
        for i in range(1, len(sand) + 1):
            if i == len(sand) or sand[i] != sand[start]:
                beds.append([float(d[start]), float(d[i - 1]) + step, bool(sand[start]), float(vsh[start:i].mean())])
                start = i
        changed = True
        while changed and len(beds) > 1:
            changed = False
            for j, b in enumerate(beds):
                if b[1] - b[0] < min_bed_ft:
                    nb = beds[j - 1] if j > 0 else beds[j + 1]
                    nb[0], nb[1] = min(nb[0], b[0]), max(nb[1], b[1])
                    del beds[j]; changed = True; break
        # re-merge same-class neighbours
        out = [beds[0]]
        for b in beds[1:]:
            if b[2] == out[-1][2]:
                out[-1][1] = b[1]
            else:
                out.append(b)
        net = sum(b[1] - b[0] for b in out if b[2])
        return {"top_ft": float(lo), "base_ft": float(hi), "gross_ft": round(float(hi - lo), 1), "net_sand_ft": round(net, 1),
                "net_to_gross": round(net / (hi - lo), 3) if hi > lo else None, "gr_clean": round(float(g_clean), 1),
                "gr_shale": round(float(g_shale), 1), "n_samples": int(k.sum()),
                "beds": [{"top_ft": round(b[0], 1), "base_ft": round(b[1], 1), "thickness_ft": round(b[1] - b[0], 1),
                          "class": "sand" if b[2] else "shale", "vsh_mean": round(b[3], 2)} for b in out]}

    res = {"method": f"linear gamma index, clean/shale lines at p{p_clean:g}/p{p_shale:g} of the gamma over the "
                     f"interval, sand where Vsh < {cutoff}, beds under {min_bed_ft} ft merged", "step_ft": round(step, 3)}
    if intervals:
        res["intervals"] = []
        for lo, hi, label in intervals:
            r = beds_for(lo, hi)
            res["intervals"].append({"label": label, **(r or {"top_ft": lo, "base_ft": hi, "note": "too few samples or flat gamma"})})
    else:
        r = beds_for(float(depth.min()), float(depth.max()))
        res["intervals"] = [{"label": "whole log", **(r or {"note": "too few samples or flat gamma"})}]
    return res


def net_sand(las_path: str, well_id: int | None = None, gr_curve: str | None = None, cutoff: float | None = None,
             min_bed_ft: float | None = None, intervals: list[list] | None = None, csv_path: str | None = None,
             overwrite: bool = False) -> dict:
    import lasio
    p = Path(las_path).expanduser()
    if not p.exists():
        raise UserError(f"{p} does not exist")
    las = lasio.read(str(p), ignore_header_errors=True)
    name = _gamma_curve(las, gr_curve)
    df = las.df()
    cfg = {**NET_SAND_DEFAULTS, **{k: v for k, v in (("cutoff", cutoff), ("min_bed_ft", min_bed_ft)) if v is not None}}
    if not 0 < cfg["cutoff"] < 1:
        raise UserError("cutoff is a Vsh fraction between 0 and 1")
    ivals = None; dwr = {}
    if intervals:
        ivals = []
        for x in intervals:
            if len(x) < 2:
                raise UserError("each interval is [top_ft, base_ft] or [top_ft, base_ft, label]")
            ivals.append((float(x[0]), float(x[1]), str(x[2]) if len(x) > 2 else f"{x[0]}-{x[1]} ft"))
    elif well_id is not None:
        cur, m, live = _gw()
        pk = cur.execute("SELECT aquifer, top_depth_ft, base_depth_ft, net_sand_ft FROM gw_picks WHERE well_id = ? "
                         "AND top_depth_ft IS NOT NULL AND base_depth_ft IS NOT NULL ORDER BY top_depth_ft", [int(well_id)]).df()
        if pk.empty:
            raise UserError(f"well {well_id} has no picks with top and base in gw_picks; give intervals")
        ivals = [(float(r.top_depth_ft), float(r.base_depth_ft), str(r.aquifer)) for r in pk.itertuples()]
        dwr = {str(r.aquifer): (None if r.net_sand_ft != r.net_sand_ft else float(r.net_sand_ft)) for r in pk.itertuples()}
    res = net_sand_from_arrays(df.index.to_numpy(float), df[name].to_numpy(float), cfg["cutoff"], cfg["min_bed_ft"],
                               cfg["p_clean"], cfg["p_shale"], ivals)
    for iv in res["intervals"]:
        if iv["label"] in dwr:
            iv["dwr_net_sand_ft"] = dwr[iv["label"]]
            if dwr[iv["label"]] is not None and iv.get("net_sand_ft") is not None:
                iv["difference_ft"] = round(iv["net_sand_ft"] - dwr[iv["label"]], 1)
    res.update({"las": str(p), "gamma_curve": name, "digitized_from_raster": "DIGITIZED FROM RASTER" in (las.other or "").upper(),
                "notes": ["A gamma cutoff is a convention: the clean and shale lines are percentiles of this log's own gamma "
                          "over the interval, so a uniformly shaly interval still reports some 'sand'. Read the beds, not just the total.",
                          "dwr_net_sand_ft is DWR's gLogThickness for the same well and aquifer: a geologist's pick from the same log, "
                          "shown for comparison, never averaged with ours.",
                          "Benchmark of these defaults against DWR's picks: docs/research/net-sand-benchmark.md (in the repo and "
                          "the release README)."]})
    if res["digitized_from_raster"]:
        res["notes"].insert(0, "This LAS was digitized from a scan: screening data. Check the curve's gate verdict in its description.")
    if csv_path:
        import csv
        out = _safe_out(csv_path, "csv", overwrite, exts=_GW_EXT)
        with out.open("w", newline="", encoding="utf-8") as f:
            w = csv.writer(f); w.writerow(["interval", "bed_top_ft", "bed_base_ft", "thickness_ft", "class", "vsh_mean"])
            for iv in res["intervals"]:
                for b in iv.get("beds", []):
                    w.writerow([iv["label"], b["top_ft"], b["base_ft"], b["thickness_ft"], b["class"], b["vsh_mean"]])
        res["csv"] = str(out)
    return res


# --- on-demand DWR (throttled, capped, from the user's machine) -------------

_DWR_LOCK = threading.Lock()
_DWR_LAST = [0.0]


def _dwr_budget_take() -> None:
    import datetime as _dt
    p = CACHE / "dwr_calls.json"
    today = _dt.date.today().isoformat()
    try:
        d = json.loads(p.read_text(encoding="utf-8"))
    except (OSError, ValueError):
        d = {}
    n = int(d.get(today, 0))
    if n >= DWR_DAILY_CAP:
        raise UserError(f"this machine has made {n} on-demand DWR requests today (cap {DWR_DAILY_CAP}); "
                        "the bulk tables still work -- try again tomorrow")
    CACHE.mkdir(parents=True, exist_ok=True)
    p.write_text(json.dumps({today: n + 1}), encoding="utf-8")


def _dwr(url: str, data: bytes | None = None, json_body: bool = False, opener=None, timeout: float = 90):
    import time
    with _DWR_LOCK:
        _dwr_budget_take()
        wait = _DWR_LAST[0] + DWR_INTERVAL_S - time.monotonic()
        if wait > 0:
            time.sleep(wait)
        _DWR_LAST[0] = time.monotonic()
    headers = dict(UA)
    if data is not None:
        headers["Content-Type"] = ("application/json; charset=utf-8" if json_body
                                   else "application/x-www-form-urlencoded")
    req = urllib.request.Request(url, data=data, headers=headers, method="POST" if data is not None else "GET")
    try:
        with (opener or urllib.request.build_opener()).open(req, timeout=timeout) as r:
            return r.read(), r.headers
    except urllib.error.HTTPError as e:
        raise UserError(f"DWR answered HTTP {e.code} for {url.split('?')[0]}")
    except (urllib.error.URLError, TimeoutError, OSError) as e:
        raise UserError(f"could not reach DWR ({e}); the bulk tables still work offline")


def normalize_denver_tool(payload: dict) -> dict:
    rows = []
    for r in payload.get("DenverLocationResult") or []:
        name = str(r.get("aquifer", "")).strip().replace("Fox Hils", "Fox Hills")
        present = str(r.get("depth_to_top", "")).strip() not in ("", "--")
        rows.append({
            "aquifer": name, "present": present,
            "depth_to_top_ft": _num_or_none(r.get("depth_to_top")),
            "depth_to_bottom_ft": _num_or_none(r.get("depth_to_bottom")),
            "elevation_top_ft": _num_or_none(r.get("elevation_top")),
            "elevation_bottom_ft": _num_or_none(r.get("elevation_bottom")),
            "net_sand_ft": _num_or_none(r.get("net_sand")),
            "annual_approp_af": _num_or_none(r.get("annual_approp")),
            "status": (str(r.get("status") or "").strip() or None) if present else None,
            "dwr_raw": r,
        })
    outs = {o.get("name"): o.get("value") for o in payload.get("outputs") or [] if isinstance(o, dict)}
    return {"aquifers": rows, "location": payload.get("location"), "basin": payload.get("basin"),
            "warning": payload.get("warning") or None,
            "county": outs.get("argOUTcounty"), "water_division": outs.get("argOUTdivision"),
            "water_district": outs.get("argOUTdistrict"),
            "aquifer_at_surface": outs.get("argOUTdenverbasin"), "dwr_tool_version": outs.get("argOUTversion"),
            "instance": payload.get("instanceData") or None}


def denver_basin(lat, lon, acres: float, elevation_ft: float | None = None) -> dict:
    import urllib.parse
    lat, lon = _check_point(lat, lon)
    acres = float(acres)
    if not 0 < acres <= 100000:
        raise UserError("acres is the land area the determination is for (e.g. 40)")
    base = f"{DWR_TOOLS}/DenverLocation/DenverLocation"
    elev_src = "given"
    if elevation_ft is None:
        body, _ = _dwr(f"{base}/GetElevationByLocation",
                       urllib.parse.urlencode({"Lat": lat, "Long": lon}).encode())
        elevation_ft = float(int(float(body.decode().strip())))
        elev_src = "DWR GetElevationByLocation (the tool's own default)"
    body, _ = _dwr(f"{base}/getLocationDataLatLong", urllib.parse.urlencode(
        {"Lattitude": lat, "Longitude": lon, "acres": acres, "elevation": elevation_ft}).encode())
    payload = json.loads(body)
    if "NOT IN COLORADO" in str(payload.get("instanceData", "")).upper():
        raise UserError("DWR's tool says the point is not in Colorado")
    out = normalize_denver_tool(payload)
    out.update({"input": {"latitude": lat, "longitude": lon, "acres": acres, "ground_elevation_ft": elevation_ft,
                          "elevation_source": elev_src},
                "source": f"{DWR_TOOLS}/DenverLocation", "note": DENVER_TOOL_NOTE})
    return out


def formation_tops(lat, lon, radius_mi: float = 2.0, dwr_tool: bool = True, acres: float = 40.0,
                   elevation_ft: float | None = None) -> dict:
    """Measured picks nearby beside DWR's tool at the point. Elevations are compared (depths
    depend on each well's own ground). Nothing is averaged or interpolated."""
    lat, lon = _check_point(lat, lon)
    cur, m, live = _gw()
    picks = cur.execute(
        f"SELECT round({_dist_sql(lat, lon, 'w.latitude', 'w.longitude')}, 3) AS distance_mi, "
        f"p.well_id, w.name, w.log_date, w.elevation_ft AS well_ground_ft, p.aquifer, p.top_depth_ft, "
        f"p.base_depth_ft, p.top_elev_ft, p.base_elev_ft, p.interval_ft, p.net_sand_ft, p.comment "
        f"FROM gw_picks p JOIN gw_log_wells w USING (well_id) "
        f"WHERE {_bbox_sql(lat, lon, float(radius_mi), 'w.latitude', 'w.longitude')}").df()
    picks = picks[picks["distance_mi"] <= float(radius_mi)].sort_values(["distance_mi", "top_depth_ft"])
    out = {"picks": _rows(picks, stale=not live),
           "pick_note": "Picks are DWR's from each well's geophysical log. net_sand_ft is DWR's 'total "
                        "thickness of the sand layers'; interval_ft is base minus top. Coarse picks "
                        "(DAWSON, ARAPAHOE) span the tool's Upper+Lower and are not compared to either."}
    if dwr_tool:
        try:
            tool = denver_basin(lat, lon, acres, elevation_ft)
        except UserError as e:
            out["dwr_tool"] = {"error": str(e)}
            return out
        out["dwr_tool"] = tool
        comp = []
        for a in tool["aquifers"]:
            near = picks[picks["aquifer"].map(lambda x: PICK_TO_TOOL.get(str(x).upper())) == a["aquifer"]]
            row = {"aquifer": a["aquifer"], "tool_top_elev_ft": a["elevation_top_ft"],
                   "tool_bottom_elev_ft": a["elevation_bottom_ft"], "tool_net_sand_ft": a["net_sand_ft"]}
            coarse = picks[picks["aquifer"].str.upper() == COARSE_PICK.get(a["aquifer"], "(none)")]
            if not coarse.empty:
                c = coarse.iloc[0]
                row["nearest_coarse_pick"] = {
                    "well_id": int(c["well_id"]), "name": c["name"], "distance_mi": float(c["distance_mi"]),
                    "aquifer": c["aquifer"], "top_elev_ft": _num_or_none(c["top_elev_ft"]),
                    "base_elev_ft": _num_or_none(c["base_elev_ft"]), "net_sand_ft": _num_or_none(c["net_sand_ft"]),
                    "note": f"picked as {c['aquifer']} (spans Upper and Lower); shown for context, not compared"}
            if near.empty:
                row["nearest_pick"] = None
            else:
                p = near.iloc[0]
                row["nearest_pick"] = {"well_id": int(p["well_id"]), "name": p["name"],
                                       "distance_mi": float(p["distance_mi"]),
                                       "top_elev_ft": _num_or_none(p["top_elev_ft"]),
                                       "base_elev_ft": _num_or_none(p["base_elev_ft"]),
                                       "net_sand_ft": _num_or_none(p["net_sand_ft"])}
                if a["elevation_top_ft"] is not None and row["nearest_pick"]["top_elev_ft"] is not None:
                    row["top_elev_difference_ft"] = round(row["nearest_pick"]["top_elev_ft"] - a["elevation_top_ft"], 1)
            comp.append(row)
        out["comparison"] = comp
        out["comparison_note"] = ("Each tool aquifer beside the NEAREST pick of the same name, by elevation. "
                                  "A difference is a disagreement between two DWR products to report, not to "
                                  "average away; the pick is a measurement at its well, the tool an "
                                  "interpolated surface at your point.")
    return out


def cross_section(well_ids: list | None = None, line: list | None = None, buffer_mi: float = 0.5,
                  csv_path: str | None = None, overwrite: bool = False) -> dict:
    """Picks along a section: either the wells you name (in that order) or every log well within
    buffer_mi of a polyline [[lat, lon], ...]. Distance is along the line; no interpolation."""
    import math
    import pandas as pd
    cur, m, live = _gw()
    if well_ids:
        ids = [int(i) for i in well_ids][:50]
        w = cur.execute(f"SELECT well_id, name, latitude, longitude, elevation_ft FROM gw_log_wells "
                        f"WHERE well_id IN ({','.join(map(str, ids))})").df().set_index("well_id").reindex(ids).dropna(
            subset=["latitude"]).reset_index()
        pts = list(zip(w["latitude"], w["longitude"]))
        along = [0.0]
        for (a, b), (c, d) in zip(pts, pts[1:]):
            along.append(along[-1] + _hav(a, b, c, d))
        w["station_mi"] = along
        w["offset_mi"] = 0.0
    elif line and len(line) >= 2:
        verts = [_check_point(a, b) for a, b in line]
        lats = [v[0] for v in verts]
        lons = [v[1] for v in verts]
        r = float(buffer_mi)
        cand = cur.execute(
            f"SELECT well_id, name, latitude, longitude, elevation_ft FROM gw_log_wells WHERE latitude BETWEEN "
            f"{min(lats) - r / 69} AND {max(lats) + r / 69} AND longitude BETWEEN {min(lons) - r / 50} AND "
            f"{max(lons) + r / 50} AND latitude IS NOT NULL").df()
        rows = []
        for t in cand.itertuples():
            st, off = _project(t.latitude, t.longitude, verts)
            if off <= r:
                rows.append({"well_id": t.well_id, "name": t.name, "latitude": t.latitude,
                             "longitude": t.longitude, "elevation_ft": t.elevation_ft,
                             "station_mi": round(st, 3), "offset_mi": round(off, 3)})
        w = pd.DataFrame(rows, columns=["well_id", "name", "latitude", "longitude", "elevation_ft",
                                        "station_mi", "offset_mi"]).sort_values("station_mi")
    else:
        raise UserError("give well_ids (in section order) or line=[[lat, lon], [lat, lon], ...]")
    if w.empty:
        raise UserError("no log wells on that section; widen buffer_mi or check gw_logs_near")
    p = cur.execute(f"SELECT well_id, aquifer, top_depth_ft, base_depth_ft, top_elev_ft, base_elev_ft, "
                    f"interval_ft, net_sand_ft FROM gw_picks WHERE well_id IN "
                    f"({','.join(str(int(i)) for i in w['well_id'])})").df()
    sec = w.merge(p, on="well_id", how="left").sort_values(["station_mi", "top_depth_ft"])
    res = {"section": _rows(sec, stale=not live),
           "note": "One row per pick per well, ordered along the section. Wells are posts, not a surface: "
                   "nothing between them is interpolated. offset_mi is each well's distance off the line."}
    if csv_path:
        out = _safe_out(csv_path, "csv", overwrite, exts=_GW_EXT)
        out.parent.mkdir(parents=True, exist_ok=True)
        sec.to_csv(out, index=False)
        res["written"] = str(out)
    return res


def _hav(a, b, c, d) -> float:
    import math
    p1, p2 = math.radians(a), math.radians(c)
    h = math.sin((p2 - p1) / 2) ** 2 + math.cos(p1) * math.cos(p2) * math.sin(math.radians(d - b) / 2) ** 2
    return 2 * EARTH_MI * math.asin(math.sqrt(h))


def _project(lat, lon, verts) -> tuple[float, float]:
    """(station along the polyline, perpendicular offset), miles, on a local flat projection."""
    import math
    k = math.cos(math.radians(lat))
    x, y = lon * 69.0 * k, lat * 69.0
    best, run = (float("inf"), 0.0), 0.0
    for (a, b), (c, d) in zip(verts, verts[1:]):
        x1, y1, x2, y2 = b * 69.0 * k, a * 69.0, d * 69.0 * k, c * 69.0
        dx, dy = x2 - x1, y2 - y1
        seg = math.hypot(dx, dy) or 1e-9
        t = max(0.0, min(1.0, ((x - x1) * dx + (y - y1) * dy) / seg ** 2))
        off = math.hypot(x - (x1 + t * dx), y - (y1 + t * dy))
        if off < best[0]:
            best = (off, run + t * seg)
        run += seg
    return best[1], best[0]


def gw_documents_list(well_id: int | None = None, receipt: str | None = None,
                      doc_type: str | None = None) -> dict:
    if well_id is None and not receipt:
        raise UserError("give well_id (a log well) or receipt")
    cur, m, live = _gw()
    where, params = [], []
    if well_id is not None:
        where.append("well_id = ?")
        params.append(int(well_id))
    if receipt:
        where.append("(receipt = ? OR well_id IN (SELECT well_id FROM gw_log_wells WHERE receipt = ?))")
        params += [str(receipt).strip()] * 2
    sql = f"SELECT * FROM gw_documents WHERE {' OR '.join(where)}"
    if doc_type:
        sql = f"SELECT * FROM ({sql}) WHERE doc_type ILIKE ?"
        params.append(f"%{doc_type}%")
    df = cur.execute(sql + " ORDER BY category, created", params).df()
    out = {"documents": _rows(df, stale=not live)}
    if df.empty:
        out["note"] = ("none listed in the release. The mirror indexes log wells and monitored wells; for "
                       "other permits, DWR's permit page (gw_permit -> dwr_url) lists every document.")
    return out


def _weblink_pdf(doc_id: int) -> bytes:
    import http.cookiejar
    import time
    op = urllib.request.build_opener(urllib.request.HTTPCookieProcessor(http.cookiejar.CookieJar()))
    _dwr(f"{WEBLINK}/DocView.aspx?id={doc_id}&dbid=0", opener=op)
    body, _ = _dwr(f"{WEBLINK}/GeneratePDF10.aspx?key={doc_id}&PageRange=&Watermark=", data=b"{}",
                   json_body=True, opener=op)
    key = body.decode("utf-8", "replace").splitlines()[0].strip()
    if not re.fullmatch(r"[0-9a-fA-F-]{36}", key):
        raise UserError(f"DWR's PDF export did not start for document {doc_id}")
    for _ in range(20):
        prog, _ = _dwr(f"{WEBLINK}/DocumentService.aspx/PDFTransition", data=json.dumps({"Key": key}).encode(),
                       json_body=True, opener=op)
        d = json.loads(prog).get("data") or {}
        if d.get("errMsg"):
            raise UserError(f"DWR's PDF export failed: {d['errMsg']}")
        if d.get("finished"):
            break
        time.sleep(1)
    else:
        raise UserError("DWR's PDF export did not finish; try again later")
    pdf, _ = _dwr(f"{WEBLINK}/PDF10/{key}/{doc_id}", opener=op, timeout=300)
    if not pdf.startswith(b"%PDF"):
        raise UserError("DWR's PDF export returned something other than a PDF")
    return pdf


def fetch_document(doc_id: int, path: str, overwrite: bool = False) -> dict:
    cur, m, live = _gw()
    df = cur.execute("SELECT * FROM gw_documents WHERE doc_id = ? LIMIT 1", [int(doc_id)]).df()
    if df.empty:
        raise UserError(f"document {doc_id} is not in the release; gw_documents lists a well's documents")
    d = df.iloc[0].to_dict()
    if d.get("withdrawn_on") and str(d["withdrawn_on"]) not in ("None", "nan", "NaT", ""):
        raise UserError(f"Hydrohistorian withdrew its copy of document {doc_id} on {d['withdrawn_on']} at a "
                        f"request (hydrohistorian.com/well-documents). DWR's own record: {d['dwr_url']}")
    is_las = str(d["doc_type"]).upper().endswith("LAS")
    want = "las" if is_las else ("pdf" if "DocView" in str(d["dwr_url"]) else (d.get("format") or "pdf"))
    out = _safe_out(path, want if want in _GW_EXT else "pdf", overwrite, exts=_GW_EXT)
    blob, via = None, None
    if d.get("mirrored") and FILES_URL and d.get("file"):
        try:
            req = urllib.request.Request(f"{FILES_URL}/{d['file']}", headers=UA)
            with urllib.request.urlopen(req, timeout=300) as r:
                blob = r.read()
            if hashlib.sha256(blob).hexdigest() != d.get("sha256"):
                blob = None
            else:
                via = "Hydrohistorian mirror (sha256 verified)"
        except Exception:  # noqa: BLE001 -- fall through to DWR
            blob = None
    if blob is None:
        if "ElectronicFile" in str(d["dwr_url"]):
            blob, _ = _dwr(f"{WEBLINK}/ElectronicFile.aspx?docid={int(doc_id)}&dbid=0", timeout=300)
        else:
            blob = _weblink_pdf(int(doc_id))
        via = "DWR Laserfiche, fetched now"
    out.parent.mkdir(parents=True, exist_ok=True)
    out.write_bytes(blob)
    return {"path": str(out), "bytes": len(blob), "doc_type": d["doc_type"], "via": via,
            "dwr_url": d["dwr_url"], "sha256": hashlib.sha256(blob).hexdigest()}


# --- LAS files (DWR's 'Log LAS', ECMC, USGS, a logging contractor's) ------------

def _las_in(path: str):
    p = Path(path).expanduser().resolve()
    if p.suffix.lower() != ".las":
        raise UserError("a LAS file ends in .las (gw_fetch_document saves DWR's that way)")
    if any(part.startswith(".") for part in p.parts[1:]):
        raise UserError("refusing to read from a hidden directory")
    if not p.exists():
        raise UserError(f"{p} does not exist")
    return p


def _read_las(p: Path):
    import lasio
    try:
        return lasio.read(str(p), ignore_header_errors=True, engine="normal")
    except Exception as e:  # noqa: BLE001
        raise UserError(f"{p.name} did not parse as LAS: {e}")


def read_las(path: str) -> dict:
    """Header and per-curve summary of one LAS file."""
    import numpy as np
    p = _las_in(path)
    las = _read_las(p)
    df = las.df()
    depth = df.index.to_numpy(dtype=float)
    well = {i.mnemonic: (str(i.value).strip(), i.unit, i.descr) for i in las.well}
    curves = []
    for c in las.curves[1:]:
        v = df[c.mnemonic].to_numpy(dtype=float) if c.mnemonic in df else np.array([])
        ok = ~np.isnan(v)
        curves.append({"mnemonic": c.mnemonic, "unit": c.unit or None, "description": c.descr or None,
                       "samples": int(ok.sum()),
                       "top_ft": float(depth[ok].min()) if ok.any() else None,
                       "bottom_ft": float(depth[ok].max()) if ok.any() else None,
                       "min": float(np.nanmin(v)) if ok.any() else None,
                       "max": float(np.nanmax(v)) if ok.any() else None,
                       "median": float(np.nanmedian(v)) if ok.any() else None})
    other = str(las.other).strip()
    return {"path": str(p), "version": str(las.version["VERS"].value) if "VERS" in las.version else None,
            "index": {"mnemonic": las.curves[0].mnemonic, "unit": las.curves[0].unit,
                      "start": well.get("STRT", (None,))[0], "stop": well.get("STOP", (None,))[0],
                      "step": well.get("STEP", (None,))[0]},
            "null_value": well.get("NULL", (None,))[0],
            "well": {k: v[0] for k, v in well.items() if v[0] and k not in ("STRT", "STOP", "STEP", "NULL")},
            "curves": curves,
            "other": other[:2000] or None,
            "digitized_from_raster": "DIGITIZED FROM RASTER" in other.upper(),
            "note": "Values as the file states them (LAS NULL -> missing). Units are the file's; nothing converted.",
            **({"nmr_k_note": "This file carries NMR-derived hydraulic conductivity (" + ", ".join(k) + "). Those are "
                "model estimates (SDR/SOE-type transforms with default constants), and the models can differ "
                "by an order of magnitude in one well; they are not aquifer-test results. Compare medians "
                "here before quoting either."} if (k := [c["mnemonic"] for c in curves
                                                        if re.search(r"HYDRAULIC|PERM|K_?(SDR|SOE|TC)", c["mnemonic"].upper())])
               else {})}


def las_to_csv(paths: list | str, out_path: str, curves: list | None = None, top_ft: float | None = None,
               bottom_ft: float | None = None, overwrite: bool = False) -> dict:
    """One CSV, indexed by depth, from one LAS or several runs of the same well. Runs are joined
    on depth (outer join, never resampled); a curve name in two runs gets a _run2 suffix."""
    import pandas as pd
    plist = [paths] if isinstance(paths, str) else list(paths)
    if not 1 <= len(plist) <= 10:
        raise UserError("give 1 to 10 LAS files")
    out = _safe_out(out_path, "csv", overwrite, exts=_GW_EXT)
    frames, units = [], {}
    for n, path in enumerate(plist, 1):
        las = _read_las(_las_in(path))
        df = las.df()
        df.index = df.index.astype(float).round(4)
        df.index.name = "depth"
        cols = [c for c in df.columns if not curves or c in curves]
        df = df[cols]
        for c in las.curves[1:]:
            if c.mnemonic in cols:
                units[c.mnemonic] = c.unit
        if frames:
            have = set().union(*[f.columns for f in frames])
            df = df.rename(columns={c: f"{c}_run{n}" for c in df.columns if c in have})
            units.update({f"{c}_run{n}": units.get(c) for c in cols if f"{c}_run{n}" in df.columns})
        frames.append(df)
    merged = pd.concat(frames, axis=1, join="outer").sort_index()
    if top_ft is not None:
        merged = merged[merged.index >= float(top_ft)]
    if bottom_ft is not None:
        merged = merged[merged.index <= float(bottom_ft)]
    merged = merged.dropna(how="all")
    if merged.empty:
        raise UserError("no samples left (check curves= names against gw_read_las, and the depth window)")
    merged.columns = [f"{c} ({units[c]})" if units.get(c) else c for c in merged.columns]
    out.parent.mkdir(parents=True, exist_ok=True)
    merged.to_csv(out)
    return {"path": str(out), "rows": int(len(merged)), "columns": ["depth"] + list(merged.columns),
            "depth_range": [float(merged.index.min()), float(merged.index.max())]}


# --- Raster log digitizer -------------------------------------------------------
#
# Semi-automatic, on purpose: the caller (you, or the model looking at the page
# with gw_log_page) supplies where depth and each track's scale are; the code
# only traces ink. Every output is labelled DIGITIZED FROM RASTER, carries its
# calibration, a per-curve coverage figure and the benchmark below, and comes
# with a QA overlay to look at before anyone quotes a value.

DIGITIZER_VERSION = "1"
NAMED_HUES = {"red": 0, "orange": 30, "yellow": 60, "green": 120, "cyan": 180, "blue": 240, "magenta": 300,
              "purple": 280}
# Measured on DWR's own scan+LAS pairs (docs/research/log-digitizer-benchmark.md); filled in by that run.
DIGITIZER_BENCHMARK = (
    "v1 vs DWR's own LAS for the same run (1 well, Colog 2026 print; "
    "https://hydrohistorian.com/static/data/releases/gw/digitizer-benchmark.md): colored curves r>=0.994, mean error 0.2-0.7% of track scale, depth lag 0; "
    "gamma r 0.84 because the print was smoothed (~1.5 ft; r 0.994 vs smoothed truth); gray photocopy, "
    "one curve per track: r 0.84-0.99, ~1% of scale; several same-ink curves in one track: unusable "
    "(flagged). Untested: skewed/warped, faded or hand-drawn scans. Screening data -- keep the scan. "
    "The same document's multi-well section scores the AUTOMATIC digitizer on 377 scan+LAS pairs; its "
    "output is the gw_digitized table (gw_digitized_curves / gw_fetch_digitized).")


def _load_page(path: str, page: int = 1):
    """RGB uint8 array of an image file, or of one page of a PDF (its embedded scan when the
    page is a single image, which is how DWR's Laserfiche exports are built)."""
    import io
    import numpy as np
    from PIL import Image
    p = Path(path).expanduser().resolve()
    if not p.exists():
        raise UserError(f"{p} does not exist")
    Image.MAX_IMAGE_PIXELS = 400_000_000  # log scans are tall (2700 x 53100 is normal)
    if p.suffix.lower() == ".pdf":
        from pypdf import PdfReader
        r = PdfReader(str(p))
        if not 1 <= int(page) <= len(r.pages):
            raise UserError(f"page is 1..{len(r.pages)}")
        imgs = list(r.pages[int(page) - 1].images)
        if not imgs:
            raise UserError("that PDF page holds no scanned image (a vector plot needs a different route)")
        im = max((Image.open(io.BytesIO(i.data)) for i in imgs), key=lambda x: x.size[0] * x.size[1])
    elif p.suffix.lower() in (".png", ".tif", ".tiff", ".jpg", ".jpeg"):
        im = Image.open(p)
    else:
        raise UserError("give a .pdf, .png, .tif or .jpg")
    return np.asarray(im.convert("RGB"))


def log_page(path: str, page: int = 1, out_png: str | None = None, y0: int = 0, y1: int | None = None,
             overwrite: bool = False) -> dict:
    """A window of a log page with a pixel ruler, so its coordinates can be read for calibration."""
    import numpy as np
    from PIL import Image, ImageDraw
    a = _load_page(path, page)
    H, W, _ = a.shape
    y1 = H if y1 is None else min(int(y1), H)
    y0 = max(0, int(y0))
    if y1 - y0 > 2400:
        raise UserError(f"window at most 2400 px tall (page is {W} x {H}); step through it")
    res = {"width": W, "height": H, "window": [y0, y1]}
    if out_png:
        out = _safe_out(out_png, "png", overwrite, exts=_GW_EXT)
        pad = 60
        canvas = Image.new("RGB", (W + pad, y1 - y0 + 20), "white")
        canvas.paste(Image.fromarray(a[y0:y1]), (pad, 20))
        d = ImageDraw.Draw(canvas)
        for y in range((y0 // 50 + 1) * 50, y1, 50):
            d.line([(pad - (12 if y % 100 == 0 else 5), y - y0 + 20), (pad, y - y0 + 20)], fill="red")
            if y % 100 == 0:
                d.text((2, y - y0 + 14), str(y), fill="red")
        for x in range(0, W, 50):
            d.line([(x + pad, 0), (x + pad, 12 if x % 100 == 0 else 6)], fill="blue")
            if x % 200 == 0:
                d.text((x + pad + 2, 2), str(x), fill="blue")
        out.parent.mkdir(parents=True, exist_ok=True)
        canvas.save(out)
        res["written"] = str(out)
        res["how_to_read"] = ("Red ticks on the left are page rows (y, every 50 px, labelled every 100); blue "
                              "ticks on top are columns (x). Give gw_digitize_log rows where you can read a depth "
                              "label on a gridline, and each track's x at its scale's left/right values.")
    return res


def _hue_sat_val(a):
    import numpy as np
    f = a.astype(np.float32) / 255.0
    mx, mn = f.max(2), f.min(2)
    d = mx - mn
    h = np.zeros_like(mx)
    r, g, b = f[..., 0], f[..., 1], f[..., 2]
    nz = d > 1e-6
    m = nz & (mx == r)
    h[m] = (60 * ((g - b)[m] / d[m])) % 360
    m = nz & (mx == g)
    h[m] = 60 * ((b - r)[m] / d[m]) + 120
    m = nz & (mx == b)
    h[m] = 60 * ((r - g)[m] / d[m]) + 240
    s = np.where(mx > 0, d / np.maximum(mx, 1e-6), 0)
    return h, s, mx


def _dark_ink(v, s):
    """Black/gray ink, thresholded per strip: halfway between the paper (median brightness) and
    the curve (the typical darkest pixel of a row). A photocopied 1970s log is gray-on-gray;
    a fixed cut-off loses its curve entirely (benchmark: SP r -0.11 -> see the doc)."""
    import numpy as np
    paper = float(np.median(v))
    ink = float(np.median(v.min(axis=1)))
    thr = ink + 0.5 * (paper - ink) if paper - ink > 0.15 else 0.35
    return (v < thr) & (s < 0.35)


def _ink_mask(strip, color: str, dark_max: float | None = None):
    """Pixels of one curve. Chromatic inks by hue (+-25 deg, saturated); 'black' by darkness.
    dark_max optionally splits two inks of one hue (dark green SPR vs bright green gamma)."""
    import numpy as np
    h, s, v = _hue_sat_val(strip)
    c = str(color).strip().lower()
    if c in ("black", "dark", "k"):
        return _dark_ink(v, s)
    if c.startswith("#") and len(c) == 7:
        rgb = np.array([[[int(c[i:i + 2], 16) for i in (1, 3, 5)]]], dtype=np.uint8)
        hh, ss, vv = _hue_sat_val(rgb)
        if ss[0, 0] < 0.25:
            return _dark_ink(v, s)
        hue = float(hh[0, 0])
    elif c in NAMED_HUES:
        hue = NAMED_HUES[c]
    else:
        raise UserError(f"color is black, a name ({', '.join(NAMED_HUES)}) or #rrggbb, not {color!r}")
    dh = np.minimum(np.abs(h - hue), 360 - np.abs(h - hue))
    m = (dh <= 25) & (s >= 0.35) & (v >= 0.15)
    if dark_max is not None:
        m &= v <= float(dark_max)
    return m


def _trace(mask, max_jump_frac: float = 0.6):
    """Column (float) of the curve in each row, NaN where no ink. Gridlines (rows or columns
    mostly 'ink') are ignored; where a row has several ink runs, the one nearest the last
    traced position wins (continuity). Returns (x, wraps) -- wraps = rows where the trace
    jumped more than max_jump_frac of the track (a backup-scale wrap or a lost curve)."""
    import numpy as np
    H, W = mask.shape
    mask = mask.copy()
    mask[:, mask.mean(0) > 0.30] = False         # vertical gridlines / track borders in this ink
    mask[mask.mean(1) > 0.60, :] = False         # horizontal gridlines in this ink
    xs = np.full(H, np.nan)
    last = None
    wraps = []
    multi = 0
    for y in range(H):
        cols = np.flatnonzero(mask[y])
        if cols.size == 0:
            continue
        breaks = np.flatnonzero(np.diff(cols) > 2)
        runs = np.split(cols, breaks + 1)
        runs = [r for r in runs if len(r) >= 2] or runs  # 1-px specks are scan noise, not a curve
        centers = np.array([(r[0] + r[-1]) / 2.0 for r in runs])
        if last is None:
            x = centers[np.argmin(np.abs(centers - np.median(cols)))]
        else:
            x = centers[np.argmin(np.abs(centers - last))]
            if abs(x - last) > max_jump_frac * W:
                wraps.append(y)
        # another real stroke of this ink well away from the curve = a second curve in the color
        multi += bool(np.any(np.abs(centers - x) > 0.05 * W))
        xs[y] = x
        last = x
    inked = int(np.sum(~np.isnan(xs)))
    return xs, wraps, (multi / inked if inked else 0.0)


def digitize_log(path: str, depth_calibration: list, tracks: list, out_las: str, page: int = 1,
                 top_row: int | None = None, bottom_row: int | None = None, step_ft: float = 0.5,
                 qa_png: str | None = None, overwrite: bool = False) -> dict:
    import numpy as np
    a = _load_page(path, page)
    H, W, _ = a.shape
    cal = np.array([[float(r), float(d)] for r, d in depth_calibration])
    if len(cal) < 2 or np.ptp(cal[:, 0]) < 100:
        raise UserError("depth_calibration needs >= 2 [row, depth] pairs at least 100 px apart "
                        "(read them off gridline labels with gw_log_page)")
    slope, icpt = np.polyfit(cal[:, 0], cal[:, 1], 1)
    resid = cal[:, 1] - (slope * cal[:, 0] + icpt)
    if len(cal) > 2 and np.abs(resid).max() > 1.0:
        raise UserError(f"the calibration points disagree by up to {np.abs(resid).max():.1f} ft from one straight "
                        "line -- check a misread label (or the scale changes on this page)")
    y0 = max(0, int(top_row or 0))
    y1 = min(H, int(bottom_row or H))
    if not tracks:
        raise UserError("give at least one track")
    out = _safe_out(out_las, "las", overwrite, exts=_GW_EXT)
    rows = np.arange(y0, y1)
    depth_rows = slope * rows + icpt
    grid = np.arange(np.ceil(depth_rows.min() / step_ft) * step_ft, depth_rows.max(), step_ft)
    curves, report, traces = {}, [], []
    for t in tracks:
        name = str(t.get("mnemonic") or t.get("name") or "").strip().upper()
        if not re.fullmatch(r"[A-Z][A-Z0-9_]{0,15}", name):
            raise UserError(f"track mnemonic {name!r}: letters/digits/underscore, <= 16")
        xl, xr = int(t["x_left"]), int(t["x_right"])
        vl, vr = float(t["value_left"]), float(t["value_right"])
        scale = str(t.get("scale", "linear")).lower()
        if not 0 <= xl < xr <= W:
            raise UserError(f"{name}: x_left/x_right must lie within 0..{W}")
        if scale == "log" and (vl <= 0 or vr <= 0):
            raise UserError(f"{name}: a log scale needs positive end values")
        strip = a[y0:y1, xl:xr + 1]
        mask = _ink_mask(strip, t.get("color", "black"), t.get("dark_max"))
        xs, wraps, ambiguous = _trace(mask)
        frac = xs / max(1, xr - xl)
        if scale == "log":
            vals = 10 ** (np.log10(vl) + frac * (np.log10(vr) - np.log10(vl)))
        else:
            vals = vl + frac * (vr - vl)
        # one value per depth step: the median of the rows that fall in it (no gap filling)
        idx = np.floor((depth_rows - grid[0]) / step_ft + 0.5).astype(int)
        out_v = np.full(len(grid), np.nan)
        ok = ~np.isnan(vals) & (idx >= 0) & (idx < len(grid))
        if ok.any():
            order = np.argsort(idx[ok], kind="stable")
            ii, vv = idx[ok][order], vals[ok][order]
            starts = np.r_[0, np.flatnonzero(np.diff(ii)) + 1]
            for s0, s1 in zip(starts, np.r_[starts[1:], len(ii)]):
                out_v[ii[s0]] = float(np.median(vv[s0:s1]))
        curves[name] = (out_v, str(t.get("unit", "")))
        cov = float(np.mean(~np.isnan(out_v))) if len(out_v) else 0.0
        rep = {"mnemonic": name, "coverage": round(cov, 3), "wrap_or_jump_rows": wraps[:20],
               "n_wraps": len(wraps), "ambiguous_rows": round(ambiguous, 3), "scale": scale, "range": [vl, vr]}
        if ambiguous > 0.10:
            rep["warning"] = (f"{ambiguous:.0%} of rows hold more than one stroke of this ink: another curve "
                              "(or heavy grid/noise) shares the color, and the trace may jump between them at "
                              "crossings. Check the QA overlay; if the curves cannot be told apart, do not use "
                              "this track.")
        report.append(rep)
        traces.append((xl, xs))
    import lasio
    las = lasio.LASFile()
    las.well["STRT"].value, las.well["STOP"].value, las.well["STEP"].value = float(grid[0]), float(grid[-1]), step_ft
    las.well["NULL"].value = -999.25
    las.append_curve("DEPT", grid, unit="FT", descr="Depth")
    for name, (v, unit) in curves.items():
        las.append_curve(name, v, unit=unit, descr=f"{name} digitized from raster")
    src = Path(path).name
    las.other = "\n".join([
        "DIGITIZED FROM RASTER -- NOT AN ORIGINAL DIGITAL LOG.",
        f"Source image: {src} page {page}; rows {y0}-{y1}.",
        f"Depth = {slope:.6f} * row + {icpt:.3f} ft, from {len(cal)} calibration points "
        f"(max residual {np.abs(resid).max():.2f} ft).",
        f"Tracks: " + "; ".join(f"{t.get('mnemonic')} x {t['x_left']}-{t['x_right']} = {t['value_left']}-"
                                  f"{t['value_right']} {t.get('scale', 'linear')} ink {t.get('color', 'black')}"
                                  for t in tracks),
        f"Each sample is the median of the traced rows within +-{step_ft / 2} ft; no gap filling.",
        f"hydrohistorian-mcp {__version__} digitizer v{DIGITIZER_VERSION}. Benchmark: {DIGITIZER_BENCHMARK}",
    ])
    out.parent.mkdir(parents=True, exist_ok=True)
    las.write(str(out), version=2.0)
    res = {"las": str(out), "depth_range_ft": [float(grid[0]), float(grid[-1])], "step_ft": step_ft,
           "ft_per_px": round(float(slope), 5), "tracks": report, "benchmark": DIGITIZER_BENCHMARK,
           "note": "Look at the QA overlay before using these values: the trace is drawn on the scan."}
    if qa_png:
        from PIL import Image, ImageDraw
        q = _safe_out(qa_png, "png", overwrite, exts=_GW_EXT)
        img = Image.fromarray(a[y0:y1].copy())
        d = ImageDraw.Draw(img)
        for xl, xs in traces:
            pts = [(xl + x, y) for y, x in enumerate(xs) if x == x]
            for p0, p1 in zip(pts, pts[1:]):
                if p1[1] - p0[1] <= 3:
                    d.line([p0, p1], fill=(255, 0, 255), width=1)
        q.parent.mkdir(parents=True, exist_ok=True)
        img.save(q)
        res["qa_png"] = str(q)
    return res


def _safe(fn, *a, **k) -> dict:
    """Every tool returns its error as data, so the model can read it and fix
    its call (mcp 2.x hides raised messages behind 'Error executing tool')."""
    try:
        return fn(*a, **k)
    except UserError as e:
        return {"error": str(e)}
    except Exception as e:  # noqa: BLE001
        return {"error": f"{type(e).__name__}: {e}"}


# ---------------------------------------------------------------------------
# The MCP surface
# ---------------------------------------------------------------------------

def build_server():
    try:  # mcp 2.x renamed FastMCP -> MCPServer (same constructor, tool() and run())
        from mcp.server.mcpserver import MCPServer as Server
    except ImportError:  # mcp 1.x
        from mcp.server.fastmcp import FastMCP as Server
    mcp = Server("hydrohistorian", instructions=(
        "Colorado water records from hydrohistorian.com/data, by water division (7 = San Juan/"
        "Dolores, incl. the San Juan's USGS gages in NM/UT/AZ and Reclamation's reservoirs). "
        "Start with hh_describe (tables, columns, caveats); hh_find for ids; hh_query for small "
        "answers (<=200 rows); hh_export for anything bigger -- it writes a file and keeps bulk "
        "data out of the conversation. A tool that fails returns {'error': ...}: read it and fix "
        "the call. " + GUIDE + " GROUNDWATER (statewide): division='gw' for hh_describe/hh_query/"
        "hh_export; the gw_* tools do the routine work -- gw_wells_near/gw_export_inventory for a site "
        "inventory, gw_hydrograph for water levels, gw_logs_near + gw_formation_tops + gw_cross_section "
        "for geology, gw_denver_basin for DWR's Denver Basin tool, gw_documents + gw_fetch_document for "
        "logs (LAS or scan) and construction/yield reports."))

    @mcp.tool()
    def hh_list_releases() -> dict:
        """Which water divisions have a release, and when each was built."""
        return _safe(list_releases)

    @mcp.tool()
    def hh_describe(division: int | str = 7) -> dict:
        """Tables, every column with its meaning, row counts, date ranges and caveats."""
        return _safe(describe, division)

    @mcp.tool()
    def hh_find(text: str, division: int | str = 7, limit: int = 25) -> dict:
        """Find gages, structures (WDIDs), reservoirs and streams; every word must match."""
        return _safe(find, text, division, limit)

    @mcp.tool()
    def hh_query(sql: str, division: int | str = 7) -> dict:
        """One read-only DuckDB SELECT over the division's tables; up to 200 rows back plus
        the full row count. Stopped after 60 s."""
        return _safe(query, sql, division)

    @mcp.tool()
    def hh_export(sql: str, path: str, format: str = "csv", division: int | str = 7,
                  units: str | None = None, overwrite: bool = False) -> dict:
        """Write a query's FULL result to a file: csv (.csv), parquet (.parquet) or riverware
        (.txt; needs columns date, value and units, one series only). Will not overwrite
        unless overwrite=true, and never writes hidden files."""
        return _safe(export, sql, path, format, division, units, overwrite)

    @mcp.tool()
    def hh_downstream(node_id: str, division: int | str = 7) -> dict:
        """Every node below a structure/gage on DWR's stream network, in order."""
        return _safe(downstream, node_id, division)

    @mcp.tool()
    def hh_documents(category: str | None = None, division: int | str = 7) -> dict:
        """Documents and datasets a RiverWare build of the basin rests on (StateMod, SJRIP,
        Reclamation, compacts, forecasts), with each link's latest check."""
        return _safe(documents, category, division)

    @mcp.tool()
    def gw_wells_near(lat: float, lon: float, radius_mi: float = 1.0, aquifer: str | None = None,
                      limit: int = 50) -> dict:
        """Every well around a point, statewide: well permits (depth, perforations, driller's static
        level and yield, dates, status), DWR monitoring wells (latest level, trend) and geophysical-log
        wells (log types, picks, LAS on file). The site well inventory a report starts with."""
        return _safe(wells_near, lat, lon, radius_mi, aquifer, limit)

    @mcp.tool()
    def gw_export_inventory(lat: float, lon: float, radius_mi: float, path: str, format: str = "xlsx",
                            aquifer: str | None = None, overwrite: bool = False) -> dict:
        """Write the gw_wells_near inventory to a report-ready file: xlsx (sheets Permits, Monitoring
        wells, Log wells, plus README provenance and a Columns dictionary), csv (one per table plus a
        README) or geojson (points, ready for GIS)."""
        return _safe(export_inventory, lat, lon, radius_mi, path, format, aquifer, overwrite)

    @mcp.tool()
    def gw_hydrograph(well_ids: list[int], csv_path: str | None = None, png_path: str | None = None,
                      overwrite: bool = False) -> dict:
        """Water-level history for 1-20 DWR monitoring wells: per-well summary and the site's trend
        (Mann-Kendall, >=20 readings over >=10 years); optionally every reading to CSV and a
        hydrograph PNG (alt text returned)."""
        return _safe(hydrograph, well_ids, csv_path, png_path, overwrite)

    @mcp.tool()
    def gw_permit(receipt: str | None = None, permit: str | None = None) -> dict:
        """One well permit by receipt or permit number (e.g. '22189-F'), with the documents listed for it."""
        return _safe(permit_fn, receipt, permit)

    @mcp.tool()
    def gw_logs_near(lat: float, lon: float, radius_mi: float = 2.0, family: str | None = None,
                     limit: int = 100) -> dict:
        """Geophysical-log wells near a point (gamma, resistivity, SP, NMR, density...), with DWR's
        formation picks and whether a digital LAS exists. family filters on decoded curve families."""
        return _safe(logs_near, lat, lon, radius_mi, family, limit)

    @mcp.tool()
    def gw_formation_tops(lat: float, lon: float, radius_mi: float = 2.0, dwr_tool: bool = True,
                          acres: float = 40.0, elevation_ft: float | None = None) -> dict:
        """Denver Basin aquifer tops/bases at a point: measured picks from nearby log wells beside
        DWR's Denver Basin tool for the point, compared by elevation; disagreements shown, never averaged."""
        return _safe(formation_tops, lat, lon, radius_mi, dwr_tool, acres, elevation_ft)

    @mcp.tool()
    def gw_denver_basin(lat: float, lon: float, acres: float, elevation_ft: float | None = None) -> dict:
        """DWR's Denver Basin location tool for a point and acreage: per aquifer depth/elevation to top
        and bottom, net sand, annual appropriation (AF) and NT/NNT status, as DWR computes them.
        Asks DWR now (one or two requests); not a State Engineer determination."""
        return _safe(denver_basin, lat, lon, acres, elevation_ft)

    @mcp.tool()
    def gw_cross_section(well_ids: list[int] | None = None, line: list[list[float]] | None = None,
                         buffer_mi: float = 0.5, csv_path: str | None = None, overwrite: bool = False) -> dict:
        """Formation picks along a section: the log wells you name in order, or every log well within
        buffer_mi of a polyline [[lat, lon], ...]. Station and offset in miles; no interpolation."""
        return _safe(cross_section, well_ids, line, buffer_mi, csv_path, overwrite)

    @mcp.tool()
    def gw_documents(well_id: int | None = None, receipt: str | None = None, doc_type: str | None = None) -> dict:
        """DWR documents for a log well or receipt: geophysical logs (scans and LAS), well construction
        and yield-test reports, pump installation, abandonment -- each with DWR's link."""
        return _safe(gw_documents_list, well_id, receipt, doc_type)

    @mcp.tool()
    def gw_fetch_document(doc_id: int, path: str, overwrite: bool = False) -> dict:
        """Save one document (doc_id from gw_documents) to disk: LAS as .las, scans as .pdf. From
        Hydrohistorian's mirror when it holds it (checksum verified), else from DWR now."""
        return _safe(fetch_document, doc_id, path, overwrite)

    @mcp.tool()
    def gw_digitized_curves(well_id: int | None = None, lat: float | None = None, lon: float | None = None,
                            radius_mi: float = 2.0, gate: str = "pass", limit: int = 200) -> dict:
        """Curves digitized automatically from the mirrored log scans (gamma, resistivity, SP...) for a
        log well or near a point, each with its confidence-gate verdict, coverage and printed scale.
        gate='pass' (default) | 'fail' (with reasons) | 'all'. Screening data: keep the scan."""
        return _safe(digitized_curves, well_id, lat, lon, radius_mi, gate, limit)

    @mcp.tool()
    def gw_fetch_digitized(doc_id: int, page: int, path: str, overwrite: bool = False) -> dict:
        """Save the auto-digitized LAS for one scan page (doc_id + page from gw_digitized_curves) to disk."""
        return _safe(fetch_digitized, doc_id, page, path, overwrite)

    @mcp.tool()
    def gw_net_sand(las_path: str, well_id: int | None = None, gr_curve: str | None = None, cutoff: float | None = None,
                    min_bed_ft: float | None = None, intervals: list[list] | None = None, csv_path: str | None = None,
                    overwrite: bool = False) -> dict:
        """Net sand from a LAS gamma curve: linear gamma index, cutoff, beds, net sand per interval.
        Intervals come from the well's DWR aquifer picks (well_id) or from you ([[top, base, label], ...]);
        DWR's own net-sand pick is shown beside ours where it exists. csv_path writes the beds."""
        return _safe(net_sand, las_path, well_id, gr_curve, cutoff, min_bed_ft, intervals, csv_path, overwrite)

    @mcp.tool()
    def gw_read_las(path: str) -> dict:
        """Header and per-curve summary (units, depth range, min/median/max) of a LAS file --
        DWR's (gw_fetch_document), ECMC's, USGS's or your logging contractor's; NMR included."""
        return _safe(read_las, path)

    @mcp.tool()
    def gw_las_to_csv(paths: list[str], out_path: str, curves: list[str] | None = None,
                      top_ft: float | None = None, bottom_ft: float | None = None,
                      overwrite: bool = False) -> dict:
        """LAS -> one depth-indexed CSV (units in the headers). Several runs of a well merge on
        depth without resampling; curves= picks mnemonics; top_ft/bottom_ft trim."""
        return _safe(las_to_csv, paths, out_path, curves, top_ft, bottom_ft, overwrite)

    @mcp.tool()
    def gw_log_page(path: str, page: int = 1, out_png: str | None = None, y0: int = 0,
                    y1: int | None = None, overwrite: bool = False) -> dict:
        """Size of a log scan page (PDF from gw_fetch_document, or PNG/TIFF), and optionally a <=2400 px
        window of it with a pixel ruler -- look at it to read depth labels and track scales for
        gw_digitize_log."""
        return _safe(log_page, path, page, out_png, y0, y1, overwrite)

    @mcp.tool()
    def gw_digitize_log(path: str, depth_calibration: list[list[float]], tracks: list[dict], out_las: str,
                        page: int = 1, top_row: int | None = None, bottom_row: int | None = None,
                        step_ft: float = 0.5, qa_png: str | None = None, overwrite: bool = False) -> dict:
        """Trace curves off a scanned log into LAS 2.0 (labelled DIGITIZED FROM RASTER).
        depth_calibration: >=2 [row, depth_ft] pairs read on labelled gridlines (far apart).
        tracks: [{mnemonic, unit, x_left, x_right, value_left, value_right, scale: linear|log,
        color: black|red|blue|green|...|#rrggbb, dark_max (optional, splits a dark from a bright ink
        of one hue)}]. Write qa_png and LOOK at it: the trace is drawn on the scan. Accuracy:
        see 'benchmark' in the result."""
        return _safe(digitize_log, path, depth_calibration, tracks, out_las, page, top_row, bottom_row,
                     step_ft, qa_png, overwrite)

    return mcp


if __name__ == "__main__":
    build_server().run()