Graham McIntire
2f0424a47b
The hrrr_climatology table previously rebuilt only when an operator remembered to run `mix hrrr_climatology` from a workstation, which the prod environment can't do. AdminTaskWorker already exposes the same logic as a "climatology" task, so a one-line cron entry (02:30 UTC daily) lets prod self-heal: empty/stale rows refold from the hrrr_profiles archive within 24h, no manual intervention. Recalibration audit remediation now points at the cron path instead of the dev-only mix task.
929 lines
36 KiB
Python
Executable file
929 lines
36 KiB
Python
Executable file
#!/usr/bin/env python3
|
||
"""
|
||
Algo recalibration analysis runner.
|
||
|
||
Connects to the prod (or any) Postgres, joins contacts to the nearest HRRR
|
||
grid profile within ±0.07° / ±1h (matching `Weather.find_nearest_hrrr/3`),
|
||
falls back to `narr_profiles` (NCEP NARR, 32 km / 3-hourly) for pre-2014
|
||
contacts where HRRR doesn't reach, and produces a Markdown report covering:
|
||
|
||
* row counts and date ranges for the key tables (HRRR + NARR + ...)
|
||
* per-band contact distribution from 50 MHz up (VHF/UHF/microwave)
|
||
* monthly sounding ducting refresh
|
||
* native HRRR profile duct stats
|
||
* per-band Pearson correlations of contact distance vs HRRR fields
|
||
* per-band Pearson correlations vs NARR fields (historical sanity check)
|
||
* per-band binned distance distributions for HPBL, pressure, NEXRAD,
|
||
and native-profile best duct band — every band with ≥50 matched
|
||
contacts gets its own bin table so we're not hard-coding the model
|
||
to 10/24 GHz behaviour
|
||
* NEXRAD composite reflectivity vs contact distance (rain-scoring check)
|
||
* hrrr_native_profiles.best_duct_band_ghz as a distance discriminator
|
||
* commercial-link rx_power degradation vs contemporaneous DFW contacts
|
||
|
||
Output is saved to docs/algo-reports/YYYY-MM-DD-recalibration.md so it can
|
||
be diffed against algo.md or committed alongside it.
|
||
|
||
Connection string is read from PROP_PROD_DB_URL (set in .envrc) or the
|
||
positional --dsn argument. Dependency-light: psycopg + pandas only.
|
||
|
||
Usage:
|
||
direnv allow # picks up PROP_PROD_DB_URL from .envrc
|
||
python3 scripts/recalibrate_algo.py
|
||
# or override:
|
||
python3 scripts/recalibrate_algo.py --dsn 'postgres://user:pw@host/db' --out path.md
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import argparse
|
||
import datetime as dt
|
||
import os
|
||
import sys
|
||
from pathlib import Path
|
||
|
||
try:
|
||
import psycopg
|
||
from psycopg.rows import dict_row
|
||
except ImportError:
|
||
sys.exit(
|
||
"psycopg is required: pip install 'psycopg[binary]>=3.1' pandas"
|
||
)
|
||
|
||
try:
|
||
import pandas as pd
|
||
except ImportError:
|
||
sys.exit("pandas is required: pip install pandas")
|
||
|
||
|
||
REPORT_HEADER = """\
|
||
# Recalibration Analysis Report — {today}
|
||
|
||
> Auto-generated by `scripts/recalibrate_algo.py`. This report is the
|
||
> input for `algo.md` updates — diff against the prior dated section to
|
||
> see what's actually moved.
|
||
|
||
Connection: `{dsn_redacted}`
|
||
Statement timeout: {stmt_timeout}
|
||
"""
|
||
|
||
|
||
# ─── SQL fragments ──────────────────────────────────────────────────────────
|
||
|
||
# (table, time_column). `time_column = None` means just count rows without
|
||
# min/max range. Tables not present in the target DB are silently skipped —
|
||
# local dev often trails prod (e.g. `propagation_scores` moved to binary files).
|
||
ROW_COUNT_SOURCES = [
|
||
("contacts", "qso_timestamp"),
|
||
("surface_observations", "observed_at"),
|
||
("soundings", "observed_at"),
|
||
("hrrr_profiles", "valid_time"),
|
||
("hrrr_native_profiles", "valid_time"),
|
||
("narr_profiles", "valid_time"),
|
||
("iemre_observations", None),
|
||
("nexrad_observations", "observed_at"),
|
||
("rtma_observations", "valid_time"),
|
||
("terrain_profiles", None),
|
||
("propagation_scores", "valid_time"),
|
||
]
|
||
|
||
CONTACTS_BY_BAND_SQL = """
|
||
SELECT band::int AS band_mhz, count(*) AS contacts,
|
||
ROUND(AVG(distance_km::numeric)) AS avg_km,
|
||
MAX(distance_km::numeric) AS max_km
|
||
FROM contacts
|
||
WHERE pos1 IS NOT NULL AND distance_km IS NOT NULL
|
||
AND distance_km < 3000 AND flagged_invalid = false AND band >= 50
|
||
GROUP BY 1 ORDER BY 1;
|
||
"""
|
||
|
||
SOUNDING_MONTHLY_SQL = """
|
||
SELECT EXTRACT(MONTH FROM observed_at)::int AS month,
|
||
count(*) AS soundings,
|
||
ROUND(100.0 * count(*) FILTER (WHERE ducting_detected) / count(*), 1) AS ducting_pct,
|
||
ROUND(AVG(min_refractivity_gradient)::numeric, 1) AS avg_min_grad,
|
||
ROUND(AVG(precipitable_water_mm)::numeric, 1) AS avg_pwat
|
||
FROM soundings
|
||
WHERE min_refractivity_gradient IS NOT NULL
|
||
GROUP BY 1 ORDER BY 1;
|
||
"""
|
||
|
||
NATIVE_DUCT_SQL = """
|
||
SELECT
|
||
CASE WHEN best_duct_band_ghz IS NULL THEN 'none'
|
||
WHEN best_duct_band_ghz < 5 THEN '<5 GHz'
|
||
WHEN best_duct_band_ghz < 15 THEN '5-15 GHz'
|
||
WHEN best_duct_band_ghz < 30 THEN '15-30 GHz'
|
||
WHEN best_duct_band_ghz < 75 THEN '30-75 GHz'
|
||
ELSE '75+ GHz' END AS band_bin,
|
||
count(*) AS profiles,
|
||
ROUND(AVG(inversion_top_m)::numeric) AS avg_inv_top_m,
|
||
ROUND(AVG(theta_e_jump_k)::numeric, 1) AS avg_theta_e_jump,
|
||
ROUND(AVG(bulk_richardson)::numeric, 2) AS avg_richardson
|
||
FROM hrrr_native_profiles GROUP BY 1 ORDER BY 1;
|
||
"""
|
||
|
||
# Per-band contact ↔ HRRR Pearson correlations (2014-10 onward).
|
||
PER_BAND_JOIN_SQL = """
|
||
WITH joined AS (
|
||
SELECT DISTINCT ON (c.id)
|
||
c.id, c.band::int AS band, c.distance_km::float AS dist,
|
||
h.surface_temp_c AS tc, h.surface_dewpoint_c AS dpc,
|
||
h.surface_pressure_mb AS pr, h.pwat_mm AS pwat,
|
||
h.min_refractivity_gradient AS grad,
|
||
h.surface_refractivity AS sref, h.hpbl_m AS hpbl
|
||
FROM contacts c
|
||
JOIN hrrr_profiles h
|
||
ON h.lat BETWEEN (c.pos1->>'lat')::float - 0.07
|
||
AND (c.pos1->>'lat')::float + 0.07
|
||
AND h.lon BETWEEN (c.pos1->>'lon')::float - 0.07
|
||
AND (c.pos1->>'lon')::float + 0.07
|
||
AND h.valid_time BETWEEN c.qso_timestamp - INTERVAL '1 hour'
|
||
AND c.qso_timestamp + INTERVAL '1 hour'
|
||
WHERE c.pos1 IS NOT NULL AND c.distance_km < 3000 AND c.flagged_invalid = false
|
||
ORDER BY c.id, ABS(EXTRACT(EPOCH FROM h.valid_time - c.qso_timestamp))
|
||
)
|
||
SELECT * FROM joined;
|
||
"""
|
||
|
||
# Per-band contact ↔ NARR Pearson correlations (pre-2014-10).
|
||
# NARR is 32 km Lambert conformal / 3-hourly, so both tolerances are wider
|
||
# than the HRRR join above. Matches `NarrClient.in_coverage?/1` — valid up
|
||
# through but excluding 2014-10-02 (HRRR takes over after that).
|
||
PER_BAND_NARR_JOIN_SQL = """
|
||
WITH joined AS (
|
||
SELECT DISTINCT ON (c.id)
|
||
c.id, c.band::int AS band, c.distance_km::float AS dist,
|
||
n.surface_temp_c AS tc, n.surface_dewpoint_c AS dpc,
|
||
n.surface_pressure_mb AS pr, n.pwat_mm AS pwat,
|
||
n.min_refractivity_gradient AS grad,
|
||
n.surface_refractivity AS sref, n.hpbl_m AS hpbl
|
||
FROM contacts c
|
||
JOIN narr_profiles n
|
||
ON n.lat BETWEEN (c.pos1->>'lat')::float - 0.25
|
||
AND (c.pos1->>'lat')::float + 0.25
|
||
AND n.lon BETWEEN (c.pos1->>'lon')::float - 0.25
|
||
AND (c.pos1->>'lon')::float + 0.25
|
||
AND n.valid_time BETWEEN c.qso_timestamp - INTERVAL '2 hours'
|
||
AND c.qso_timestamp + INTERVAL '2 hours'
|
||
WHERE c.pos1 IS NOT NULL AND c.distance_km < 3000
|
||
AND c.flagged_invalid = false
|
||
AND c.qso_timestamp < TIMESTAMP '2014-10-02'
|
||
ORDER BY c.id, ABS(EXTRACT(EPOCH FROM n.valid_time - c.qso_timestamp))
|
||
)
|
||
SELECT * FROM joined;
|
||
"""
|
||
|
||
# NEXRAD ↔ contact join. Very tight spatial tolerance because nexrad_observations
|
||
# is a per-contact enrichment table, not a grid sample.
|
||
NEXRAD_JOIN_SQL = """
|
||
WITH joined AS (
|
||
SELECT DISTINCT ON (c.id)
|
||
c.id, c.band::int AS band, c.distance_km::float AS dist,
|
||
n.max_reflectivity_dbz AS max_dbz,
|
||
n.mean_reflectivity_dbz AS mean_dbz
|
||
FROM contacts c
|
||
JOIN nexrad_observations n
|
||
ON n.lat BETWEEN (c.pos1->>'lat')::float - 0.1
|
||
AND (c.pos1->>'lat')::float + 0.1
|
||
AND n.lon BETWEEN (c.pos1->>'lon')::float - 0.1
|
||
AND (c.pos1->>'lon')::float + 0.1
|
||
AND n.observed_at BETWEEN c.qso_timestamp - INTERVAL '15 min'
|
||
AND c.qso_timestamp + INTERVAL '15 min'
|
||
WHERE c.pos1 IS NOT NULL AND c.distance_km < 3000 AND c.flagged_invalid = false
|
||
ORDER BY c.id, ABS(EXTRACT(EPOCH FROM n.observed_at - c.qso_timestamp))
|
||
)
|
||
SELECT * FROM joined;
|
||
"""
|
||
|
||
# Native HRRR profile ↔ contact join — gives us best_duct_band_ghz per contact.
|
||
NATIVE_DUCT_JOIN_SQL = """
|
||
WITH joined AS (
|
||
SELECT DISTINCT ON (c.id)
|
||
c.id, c.band::int AS band, c.distance_km::float AS dist,
|
||
n.best_duct_band_ghz AS duct_ghz,
|
||
n.inversion_top_m AS inv_top_m,
|
||
n.bulk_richardson AS richardson
|
||
FROM contacts c
|
||
JOIN hrrr_native_profiles n
|
||
ON n.lat BETWEEN (c.pos1->>'lat')::float - 0.07
|
||
AND (c.pos1->>'lat')::float + 0.07
|
||
AND n.lon BETWEEN (c.pos1->>'lon')::float - 0.07
|
||
AND (c.pos1->>'lon')::float + 0.07
|
||
AND n.valid_time BETWEEN c.qso_timestamp - INTERVAL '1 hour'
|
||
AND c.qso_timestamp + INTERVAL '1 hour'
|
||
WHERE c.pos1 IS NOT NULL AND c.distance_km < 3000 AND c.flagged_invalid = false
|
||
ORDER BY c.id, ABS(EXTRACT(EPOCH FROM n.valid_time - c.qso_timestamp))
|
||
)
|
||
SELECT * FROM joined;
|
||
"""
|
||
|
||
# Commercial-link degradation ↔ contact join. For each contact in the DFW zone
|
||
# (within ~150 km of the Princeton TX cluster at 33.2N, -96.5W), find the
|
||
# nearest commercial sample within ±15 min and compute the dB delta from the
|
||
# 7-day trailing baseline. The boost factor in Scorer.commercial_link_boost/2
|
||
# expects the same quantity, so this query answers "does our commercial-link
|
||
# signal actually correlate with longer QSO distances when we have both?"
|
||
DFW_CLUSTER_LAT = 33.2
|
||
DFW_CLUSTER_LON = -96.5
|
||
DFW_RADIUS_DEG = 1.5 # roughly 150 km at 33° latitude
|
||
|
||
COMMERCIAL_JOIN_SQL = f"""
|
||
WITH contact_zone AS (
|
||
SELECT c.id, c.qso_timestamp, c.band::int AS band, c.distance_km::float AS dist
|
||
FROM contacts c
|
||
WHERE c.pos1 IS NOT NULL AND c.distance_km < 3000 AND c.flagged_invalid = false
|
||
AND (c.pos1->>'lat')::float BETWEEN {DFW_CLUSTER_LAT - DFW_RADIUS_DEG}
|
||
AND {DFW_CLUSTER_LAT + DFW_RADIUS_DEG}
|
||
AND (c.pos1->>'lon')::float BETWEEN {DFW_CLUSTER_LON - DFW_RADIUS_DEG}
|
||
AND {DFW_CLUSTER_LON + DFW_RADIUS_DEG}
|
||
),
|
||
current_sample AS (
|
||
SELECT DISTINCT ON (cz.id)
|
||
cz.id, cz.qso_timestamp, cz.band, cz.dist,
|
||
cs.link_id, cs.rx_power_0 AS current_rx, cs.sampled_at
|
||
FROM contact_zone cz
|
||
JOIN commercial_samples cs
|
||
ON cs.sampled_at BETWEEN cz.qso_timestamp - INTERVAL '15 min'
|
||
AND cz.qso_timestamp + INTERVAL '15 min'
|
||
AND cs.link_state = 1 AND cs.rx_power_0 IS NOT NULL
|
||
ORDER BY cz.id, ABS(EXTRACT(EPOCH FROM cs.sampled_at - cz.qso_timestamp))
|
||
),
|
||
baseline AS (
|
||
SELECT cs.id, cs.qso_timestamp, cs.band, cs.dist, cs.link_id, cs.current_rx,
|
||
(
|
||
SELECT AVG(bs.rx_power_0)
|
||
FROM commercial_samples bs
|
||
WHERE bs.link_id = cs.link_id
|
||
AND bs.link_state = 1
|
||
AND bs.rx_power_0 IS NOT NULL
|
||
AND bs.sampled_at BETWEEN cs.sampled_at - INTERVAL '7 days'
|
||
AND cs.sampled_at - INTERVAL '15 min'
|
||
) AS baseline_rx
|
||
FROM current_sample cs
|
||
)
|
||
SELECT id, band, dist,
|
||
(baseline_rx - current_rx)::float AS degradation_db,
|
||
current_rx::float AS current_dbm,
|
||
baseline_rx::float AS baseline_dbm
|
||
FROM baseline
|
||
WHERE baseline_rx IS NOT NULL;
|
||
"""
|
||
|
||
# Empty-table / data-quality smoke checks. Each metric is classified in
|
||
# DATA_GAP_RULES below so the report can flag BROKEN vs WARN vs OK with
|
||
# concrete remediation hints instead of dumping a flat key/value list.
|
||
DATA_GAP_SQL = """
|
||
SELECT
|
||
(SELECT count(*) FROM narr_profiles) AS narr_rows,
|
||
(SELECT count(*) FROM narr_profiles WHERE valid_time < '2014-10-02') AS narr_pre2014_rows,
|
||
(SELECT count(*) FROM hrrr_climatology) AS hrrr_climatology_rows,
|
||
(SELECT count(*) FROM rtma_observations) AS rtma_rows,
|
||
(SELECT count(*) FROM metar_5min_observations) AS metar5_rows,
|
||
(SELECT count(*) FROM contacts WHERE hrrr_status = 'complete') AS hrrr_complete_contacts,
|
||
(SELECT count(*) FROM contacts WHERE hrrr_status != 'complete') AS hrrr_pending_contacts,
|
||
(SELECT count(*) FROM contacts
|
||
WHERE qso_timestamp < '2014-10-02'
|
||
AND pos1 IS NOT NULL AND distance_km < 3000
|
||
AND flagged_invalid = false) AS pre2014_contacts;
|
||
"""
|
||
|
||
# Coverage breakdown for tables whose total row count alone hides the
|
||
# real gap — e.g. narr_profiles with 358 rows could be one good year or
|
||
# 358 stragglers spread across 30. Year-bucketed counts let the report
|
||
# flag whether the backfill ever made meaningful progress.
|
||
NARR_COVERAGE_BY_YEAR_SQL = """
|
||
SELECT EXTRACT(YEAR FROM valid_time)::int AS year, count(*)::bigint AS profiles
|
||
FROM narr_profiles
|
||
GROUP BY 1 ORDER BY 1;
|
||
"""
|
||
|
||
# Where is the HRRR enrichment backlog? A 36k pending count could mean
|
||
# "the backfill stalled in 2018" or "every contact submitted in the last
|
||
# week is queued". The breakdown drives a different remediation: rerun
|
||
# HrrrFetchWorker.enqueue_for_qso/1 over a year range vs. unstick the
|
||
# Oban queue.
|
||
HRRR_STATUS_BY_YEAR_SQL = """
|
||
SELECT EXTRACT(YEAR FROM qso_timestamp)::int AS year,
|
||
count(*) FILTER (WHERE hrrr_status = 'complete')::bigint AS complete,
|
||
count(*) FILTER (WHERE hrrr_status = 'queued')::bigint AS queued,
|
||
count(*) FILTER (WHERE hrrr_status = 'failed')::bigint AS failed,
|
||
count(*) FILTER (WHERE hrrr_status IS NULL)::bigint AS unset
|
||
FROM contacts
|
||
WHERE pos1 IS NOT NULL AND distance_km < 3000 AND flagged_invalid = false
|
||
AND qso_timestamp >= '2014-10-02'
|
||
GROUP BY 1 ORDER BY 1;
|
||
"""
|
||
|
||
# Severity rules for each metric in DATA_GAP_SQL.
|
||
# key metric name returned by DATA_GAP_SQL
|
||
# broken_when lambda(value) → True means BROKEN (red)
|
||
# warn_when lambda(value) → True means WARN (yellow). Skipped if BROKEN.
|
||
# remediation one-line action the operator should take
|
||
# A metric not flagged by either lambda is OK (green).
|
||
DATA_GAP_RULES: list[dict] = [
|
||
{
|
||
"key": "narr_rows",
|
||
"broken_when": lambda v: v == 0,
|
||
# 358 rows is the symptom that triggered this rewrite — anything
|
||
# under ~5k means the NARR backfill hasn't meaningfully started.
|
||
"warn_when": lambda v: v < 5_000,
|
||
"remediation": "Run `mix narr.backfill` (see lib/microwaveprop/workers/narr_fetch_worker.ex).",
|
||
},
|
||
{
|
||
"key": "narr_pre2014_rows",
|
||
"broken_when": lambda v: v == 0,
|
||
"warn_when": lambda v: v < 5_000,
|
||
"remediation": (
|
||
"NARR's purpose is filling the pre-2014-10 gap; without rows here "
|
||
"the historical-sanity correlation table is empty."
|
||
),
|
||
},
|
||
{
|
||
"key": "hrrr_climatology_rows",
|
||
"broken_when": lambda v: v == 0,
|
||
"warn_when": lambda _v: False,
|
||
"remediation": (
|
||
"Should self-heal within 24h via the daily 02:30Z cron that "
|
||
"runs AdminTaskWorker(task=climatology). If still 0 after a "
|
||
"day, check Oban dashboard for failed admin jobs."
|
||
),
|
||
},
|
||
{
|
||
"key": "rtma_rows",
|
||
"broken_when": lambda v: v == 0,
|
||
"warn_when": lambda _v: False,
|
||
"remediation": (
|
||
"Enable the rtma queue and let RtmaFetchWorker run — needed for "
|
||
"contacts whose HRRR enrichment misses the ±1 h window."
|
||
),
|
||
},
|
||
{
|
||
"key": "metar5_rows",
|
||
"broken_when": lambda v: v == 0,
|
||
"warn_when": lambda _v: False,
|
||
"remediation": (
|
||
"No 5-minute METAR ingestor exists yet; if this stays at 0 the "
|
||
"schema metar_5min_observations is unused — consider dropping it "
|
||
"or wiring an IEM 5-min ASOS fetcher."
|
||
),
|
||
},
|
||
{
|
||
"key": "hrrr_complete_contacts",
|
||
# Healthy is the fraction complete, not the absolute number — handled
|
||
# separately via _check_hrrr_completeness below.
|
||
"broken_when": lambda _v: False,
|
||
"warn_when": lambda _v: False,
|
||
"remediation": "",
|
||
},
|
||
{
|
||
"key": "hrrr_pending_contacts",
|
||
"broken_when": lambda _v: False,
|
||
# Anything over 5k pending is a real backlog worth chasing.
|
||
"warn_when": lambda v: v > 5_000,
|
||
"remediation": (
|
||
"See `HRRR enrichment status by year` table below to locate the "
|
||
"backlog, then `mix hrrr.backfill --year YYYY` (or rerun "
|
||
"HrrrFetchWorker.enqueue_for_qso/1 over the affected ids)."
|
||
),
|
||
},
|
||
{
|
||
"key": "pre2014_contacts",
|
||
# pre-2014 contacts existing is informational, not a gap on its own.
|
||
"broken_when": lambda _v: False,
|
||
"warn_when": lambda _v: False,
|
||
"remediation": "",
|
||
},
|
||
]
|
||
|
||
|
||
# ─── Reporting helpers ──────────────────────────────────────────────────────
|
||
|
||
|
||
def fetch_df(conn, sql: str) -> pd.DataFrame:
|
||
with conn.cursor(row_factory=dict_row) as cur:
|
||
cur.execute(sql)
|
||
rows = cur.fetchall()
|
||
return pd.DataFrame(rows)
|
||
|
||
|
||
def md_table(df: pd.DataFrame) -> str:
|
||
if df.empty:
|
||
return "_(no rows)_\n"
|
||
return df.to_markdown(index=False, floatfmt=".3f") + "\n"
|
||
|
||
|
||
def fetch_row_counts(conn) -> pd.DataFrame:
|
||
"""Per-table counts + date ranges, tolerant of missing tables."""
|
||
rows = []
|
||
for tbl, time_col in ROW_COUNT_SOURCES:
|
||
with conn.cursor(row_factory=dict_row) as cur:
|
||
cur.execute("SELECT to_regclass(%s) AS oid", (f"public.{tbl}",))
|
||
present = cur.fetchone()["oid"] is not None
|
||
if not present:
|
||
rows.append({"tbl": tbl, "n": None, "lo": None, "hi": None})
|
||
continue
|
||
if time_col:
|
||
sql = (
|
||
f"SELECT '{tbl}' tbl, count(*)::bigint n, "
|
||
f"min({time_col})::date lo, max({time_col})::date hi FROM {tbl}"
|
||
)
|
||
else:
|
||
sql = (
|
||
f"SELECT '{tbl}' tbl, count(*)::bigint n, "
|
||
f"NULL::date lo, NULL::date hi FROM {tbl}"
|
||
)
|
||
with conn.cursor(row_factory=dict_row) as cur:
|
||
cur.execute(sql)
|
||
rows.append(cur.fetchone())
|
||
return pd.DataFrame(rows).sort_values("tbl")
|
||
|
||
|
||
def band_label(band_mhz: int) -> str:
|
||
"""Human-readable band label, matching `BandConfig` conventions."""
|
||
if band_mhz < 1_000:
|
||
return f"{band_mhz} MHz"
|
||
ghz = band_mhz / 1_000
|
||
if ghz == int(ghz):
|
||
return f"{int(ghz)} GHz"
|
||
return f"{ghz:g} GHz"
|
||
|
||
|
||
def bands_with_samples(df: pd.DataFrame, min_samples: int = 50) -> list[int]:
|
||
"""Every band in `df` with at least `min_samples` rows, sorted ascending."""
|
||
if df.empty or "band" not in df.columns:
|
||
return []
|
||
counts = df.groupby("band").size()
|
||
return sorted(int(b) for b, n in counts.items() if n >= min_samples)
|
||
|
||
|
||
def correlations_per_band(df: pd.DataFrame, min_samples: int = 30) -> pd.DataFrame:
|
||
fields = ["pr", "dpc", "pwat", "sref", "grad", "tc", "hpbl"]
|
||
rows = []
|
||
for band, group in df.groupby("band"):
|
||
if len(group) < min_samples:
|
||
continue
|
||
row = {"band_mhz": int(band), "n": len(group)}
|
||
for f in fields:
|
||
try:
|
||
row[f"rho_{f}"] = round(group["dist"].corr(group[f]), 3)
|
||
except Exception:
|
||
row[f"rho_{f}"] = None
|
||
rows.append(row)
|
||
return pd.DataFrame(rows).sort_values("band_mhz")
|
||
|
||
|
||
def hpbl_bins(df: pd.DataFrame, band_mhz: int = 10000) -> pd.DataFrame:
|
||
sub = df[(df["band"] == band_mhz) & df["hpbl"].notna()].copy()
|
||
if sub.empty:
|
||
return pd.DataFrame()
|
||
bins = [0, 200, 500, 1000, 1500, 2000, 1e9]
|
||
labels = ["<200", "200-500", "500-1000", "1000-1500", "1500-2000", ">=2000"]
|
||
sub["bin"] = pd.cut(sub["hpbl"], bins=bins, labels=labels, right=False)
|
||
out = (
|
||
sub.groupby("bin", observed=True)["dist"]
|
||
.agg(["count", "mean", "median"])
|
||
.round(1)
|
||
.rename(columns={"count": "n", "mean": "avg_km", "median": "p50_km"})
|
||
.reset_index()
|
||
)
|
||
return out
|
||
|
||
|
||
def pressure_bins(df: pd.DataFrame, band_mhz: int = 10000) -> pd.DataFrame:
|
||
sub = df[(df["band"] == band_mhz) & df["pr"].notna()].copy()
|
||
if sub.empty:
|
||
return pd.DataFrame()
|
||
bins = [0, 990, 1000, 1010, 1020, 9999]
|
||
labels = ["<990", "990-1000", "1000-1010", "1010-1020", ">=1020"]
|
||
sub["bin"] = pd.cut(sub["pr"], bins=bins, labels=labels, right=False)
|
||
out = (
|
||
sub.groupby("bin", observed=True)["dist"]
|
||
.agg(["count", "mean", "median"])
|
||
.round(1)
|
||
.rename(columns={"count": "n", "mean": "avg_km", "median": "p50_km"})
|
||
.reset_index()
|
||
)
|
||
return out
|
||
|
||
|
||
def nexrad_bins(df: pd.DataFrame, band_mhz: int = 10000) -> pd.DataFrame:
|
||
"""Distance distribution by NEXRAD max reflectivity bin, per band.
|
||
|
||
Confirms the rain-attenuation direction of the scoring: as max_dbz rises,
|
||
distance should fall at rain-sensitive bands (24+ GHz). Below-rain dBZ
|
||
(<5) is the no-precip baseline.
|
||
"""
|
||
sub = df[(df["band"] == band_mhz) & df["max_dbz"].notna()].copy()
|
||
if sub.empty:
|
||
return pd.DataFrame()
|
||
bins = [-1e9, 5, 20, 30, 45, 1e9]
|
||
labels = ["none (<5)", "drizzle (5-20)", "light (20-30)", "moderate (30-45)", "heavy (>=45)"]
|
||
sub["bin"] = pd.cut(sub["max_dbz"], bins=bins, labels=labels, right=False)
|
||
out = (
|
||
sub.groupby("bin", observed=True)["dist"]
|
||
.agg(["count", "mean", "median"])
|
||
.round(1)
|
||
.rename(columns={"count": "n", "mean": "avg_km", "median": "p50_km"})
|
||
.reset_index()
|
||
)
|
||
return out
|
||
|
||
|
||
def native_duct_bins(df: pd.DataFrame, band_mhz: int = 10000) -> pd.DataFrame:
|
||
"""Distance distribution by native-profile best duct band, per band.
|
||
|
||
Checks whether hrrr_native_profiles.best_duct_band_ghz actually predicts
|
||
longer QSO distances at the target band — the assumption behind
|
||
Scorer.score_refractivity/4's 1.15× native-duct boost.
|
||
"""
|
||
sub = df[df["band"] == band_mhz].copy()
|
||
if sub.empty:
|
||
return pd.DataFrame()
|
||
|
||
def label(v: object) -> str:
|
||
if pd.isna(v):
|
||
return "none"
|
||
v = float(v)
|
||
if v < 5:
|
||
return "<5 GHz"
|
||
if v < 10:
|
||
return "5-10 GHz"
|
||
if v < 24:
|
||
return "10-24 GHz"
|
||
if v < 47:
|
||
return "24-47 GHz"
|
||
return ">=47 GHz"
|
||
|
||
sub["bin"] = sub["duct_ghz"].map(label)
|
||
order = ["none", "<5 GHz", "5-10 GHz", "10-24 GHz", "24-47 GHz", ">=47 GHz"]
|
||
out = (
|
||
sub.groupby("bin")["dist"]
|
||
.agg(["count", "mean", "median"])
|
||
.round(1)
|
||
.rename(columns={"count": "n", "mean": "avg_km", "median": "p50_km"})
|
||
.reindex(order)
|
||
.dropna(subset=["n"])
|
||
.reset_index()
|
||
)
|
||
return out
|
||
|
||
|
||
def commercial_degradation_summary(df: pd.DataFrame) -> pd.DataFrame:
|
||
"""Does commercial-link rx_power degradation correlate with longer contacts?
|
||
|
||
Returns a per-band bin table. Empty today because commercial_samples only
|
||
covers ~2 weeks and contacts only cluster in Aug/Sep contests — but the
|
||
query is here so the moment contest-season samples land the signal shows up.
|
||
"""
|
||
if df.empty:
|
||
return pd.DataFrame()
|
||
bins = [-1e9, 0, 3, 8, 1e9]
|
||
labels = ["negative/0 (no fade)", "3 dB noise floor", "3-8 dB mild", ">=8 dB strong"]
|
||
df = df.copy()
|
||
df["bin"] = pd.cut(df["degradation_db"], bins=bins, labels=labels, right=False)
|
||
out = (
|
||
df.groupby("bin", observed=True)["dist"]
|
||
.agg(["count", "mean", "median"])
|
||
.round(1)
|
||
.rename(columns={"count": "n", "mean": "avg_km", "median": "p50_km"})
|
||
.reset_index()
|
||
)
|
||
return out
|
||
|
||
|
||
def classify_data_gaps(gaps: dict) -> pd.DataFrame:
|
||
"""Apply DATA_GAP_RULES to the raw gap counts.
|
||
|
||
Returns a dataframe with columns: metric, value, severity, remediation.
|
||
`severity` is one of "BROKEN", "WARN", "OK" (informational metrics with
|
||
no rule are also "OK"). The dataframe order matches DATA_GAP_RULES so
|
||
related metrics group naturally in the rendered table.
|
||
"""
|
||
rows = []
|
||
for rule in DATA_GAP_RULES:
|
||
key = rule["key"]
|
||
if key not in gaps:
|
||
continue
|
||
value = gaps[key]
|
||
if rule["broken_when"](value):
|
||
severity = "BROKEN"
|
||
elif rule["warn_when"](value):
|
||
severity = "WARN"
|
||
else:
|
||
severity = "OK"
|
||
rows.append(
|
||
{
|
||
"metric": key,
|
||
"value": int(value) if value is not None else None,
|
||
"severity": severity,
|
||
"remediation": rule["remediation"] if severity != "OK" else "",
|
||
}
|
||
)
|
||
# Append any metric DATA_GAP_SQL returned that we don't have a rule for,
|
||
# so the report stays exhaustive even if SQL changes ahead of the rules.
|
||
handled = {r["key"] for r in DATA_GAP_RULES}
|
||
for key, value in gaps.items():
|
||
if key in handled:
|
||
continue
|
||
rows.append(
|
||
{
|
||
"metric": key,
|
||
"value": int(value) if value is not None else None,
|
||
"severity": "OK",
|
||
"remediation": "",
|
||
}
|
||
)
|
||
return pd.DataFrame(rows)
|
||
|
||
|
||
def hrrr_completeness_pct(gaps: dict) -> float | None:
|
||
"""Fraction of eligible contacts with hrrr_status = 'complete'."""
|
||
complete = gaps.get("hrrr_complete_contacts") or 0
|
||
pending = gaps.get("hrrr_pending_contacts") or 0
|
||
total = complete + pending
|
||
if total == 0:
|
||
return None
|
||
return round(100.0 * complete / total, 1)
|
||
|
||
|
||
def commercial_correlation(df: pd.DataFrame) -> float | None:
|
||
"""Pearson correlation of degradation_db with contact distance. None if n<5."""
|
||
if df.empty or len(df) < 5:
|
||
return None
|
||
return round(df["dist"].corr(df["degradation_db"]), 3)
|
||
|
||
|
||
def redact(dsn: str) -> str:
|
||
if "@" not in dsn:
|
||
return dsn
|
||
head, tail = dsn.split("@", 1)
|
||
if "://" in head and ":" in head.split("://", 1)[1]:
|
||
scheme, rest = head.split("://", 1)
|
||
user = rest.split(":", 1)[0]
|
||
return f"{scheme}://{user}:***@{tail}"
|
||
return dsn
|
||
|
||
|
||
# ─── Main ───────────────────────────────────────────────────────────────────
|
||
|
||
|
||
def main() -> int:
|
||
parser = argparse.ArgumentParser(description=__doc__)
|
||
parser.add_argument(
|
||
"--dsn",
|
||
default=os.environ.get("PROP_PROD_DB_URL"),
|
||
help="Postgres connection string (defaults to $PROP_PROD_DB_URL)",
|
||
)
|
||
parser.add_argument(
|
||
"--out",
|
||
default=None,
|
||
help="Output Markdown path (defaults to docs/algo-reports/YYYY-MM-DD-recalibration.md)",
|
||
)
|
||
parser.add_argument(
|
||
"--statement-timeout",
|
||
default="20min",
|
||
help="Postgres statement_timeout (default: 20min)",
|
||
)
|
||
args = parser.parse_args()
|
||
|
||
if not args.dsn:
|
||
print(
|
||
"error: PROP_PROD_DB_URL not set and --dsn not given; populate "
|
||
"via .envrc",
|
||
file=sys.stderr,
|
||
)
|
||
return 2
|
||
|
||
today = dt.date.today().isoformat()
|
||
out_path = Path(
|
||
args.out
|
||
or Path(__file__).resolve().parent.parent
|
||
/ "docs"
|
||
/ "algo-reports"
|
||
/ f"{today}-recalibration.md"
|
||
)
|
||
out_path.parent.mkdir(parents=True, exist_ok=True)
|
||
|
||
print(f"connecting to {redact(args.dsn)}", file=sys.stderr)
|
||
|
||
sections: list[str] = []
|
||
sections.append(
|
||
REPORT_HEADER.format(
|
||
today=today,
|
||
dsn_redacted=redact(args.dsn),
|
||
stmt_timeout=args.statement_timeout,
|
||
)
|
||
)
|
||
|
||
with psycopg.connect(args.dsn, autocommit=True) as conn:
|
||
with conn.cursor() as cur:
|
||
cur.execute(f"SET statement_timeout = '{args.statement_timeout}';")
|
||
|
||
print("• row counts", file=sys.stderr)
|
||
sections.append("\n## Row counts and date ranges\n")
|
||
sections.append(md_table(fetch_row_counts(conn)))
|
||
|
||
print("• data gaps", file=sys.stderr)
|
||
gaps = fetch_df(conn, DATA_GAP_SQL).iloc[0].to_dict()
|
||
gap_table = classify_data_gaps(gaps)
|
||
completeness = hrrr_completeness_pct(gaps)
|
||
|
||
broken = gap_table[gap_table["severity"] == "BROKEN"]
|
||
for _, row in broken.iterrows():
|
||
print(
|
||
f" ⚠ BROKEN: {row['metric']}={row['value']} — {row['remediation']}",
|
||
file=sys.stderr,
|
||
)
|
||
|
||
sections.append("\n## Data gaps and ingestion health\n")
|
||
sections.append(
|
||
"Severity legend: **BROKEN** = table is empty and feeds an "
|
||
"active scoring path; **WARN** = populated but coverage is "
|
||
"low enough to skew per-band statistics; **OK** = within "
|
||
"spec or informational. Each non-OK row carries a one-line "
|
||
"remediation hint.\n"
|
||
)
|
||
sections.append(md_table(gap_table))
|
||
if completeness is not None:
|
||
sections.append(
|
||
f"\nHRRR enrichment completeness: **{completeness}%** "
|
||
f"({gaps.get('hrrr_complete_contacts', 0):,} complete / "
|
||
f"{gaps.get('hrrr_pending_contacts', 0):,} pending).\n"
|
||
)
|
||
|
||
print("• narr coverage by year", file=sys.stderr)
|
||
narr_year_df = fetch_df(conn, NARR_COVERAGE_BY_YEAR_SQL)
|
||
sections.append("\n### NARR coverage by year\n")
|
||
if narr_year_df.empty:
|
||
sections.append(
|
||
"_narr_profiles is empty. Run `mix narr.backfill` to "
|
||
"populate the pre-2014 corpus._\n"
|
||
)
|
||
else:
|
||
sections.append(
|
||
"Year-bucketed counts. NARR is supposed to span 1979 → "
|
||
"2014-10 at 3-hourly cadence (~2,920 timestamps per year "
|
||
"per grid cell). Years far below the others are partial "
|
||
"backfills the worker hasn't finished.\n"
|
||
)
|
||
sections.append(md_table(narr_year_df))
|
||
|
||
print("• hrrr status by year", file=sys.stderr)
|
||
hrrr_year_df = fetch_df(conn, HRRR_STATUS_BY_YEAR_SQL)
|
||
sections.append("\n### HRRR enrichment status by year\n")
|
||
if hrrr_year_df.empty:
|
||
sections.append("_(no contacts with `pos1` set after 2014-10-02)_\n")
|
||
else:
|
||
sections.append(
|
||
"Locate the pending backlog: a year with high `queued`/"
|
||
"`failed`/`unset` counts is where HrrrFetchWorker needs "
|
||
"to be rerun. Recent years dominating `queued` usually "
|
||
"means the live Oban queue is stuck rather than a "
|
||
"historical backfill gap.\n"
|
||
)
|
||
sections.append(md_table(hrrr_year_df))
|
||
|
||
print("• contacts by band", file=sys.stderr)
|
||
sections.append("\n## Contacts by band (≥50 MHz)\n")
|
||
sections.append(md_table(fetch_df(conn, CONTACTS_BY_BAND_SQL)))
|
||
|
||
print("• sounding monthly", file=sys.stderr)
|
||
sections.append("\n## Monthly sounding ducting\n")
|
||
sections.append(md_table(fetch_df(conn, SOUNDING_MONTHLY_SQL)))
|
||
|
||
print("• native HRRR ducts", file=sys.stderr)
|
||
sections.append(
|
||
"\n## HRRR native-profile duct distribution (best supportable band)\n"
|
||
)
|
||
sections.append(md_table(fetch_df(conn, NATIVE_DUCT_SQL)))
|
||
|
||
print("• per-band correlations — HRRR (this is the slow one)", file=sys.stderr)
|
||
joined = fetch_df(conn, PER_BAND_JOIN_SQL)
|
||
|
||
print("• per-band correlations — NARR (pre-2014)", file=sys.stderr)
|
||
joined_narr = fetch_df(conn, PER_BAND_NARR_JOIN_SQL)
|
||
|
||
print("• NEXRAD ↔ contacts", file=sys.stderr)
|
||
nexrad = fetch_df(conn, NEXRAD_JOIN_SQL)
|
||
|
||
print("• native-duct ↔ contacts", file=sys.stderr)
|
||
native = fetch_df(conn, NATIVE_DUCT_JOIN_SQL)
|
||
|
||
print("• commercial-link ↔ contacts (DFW zone only)", file=sys.stderr)
|
||
commercial = fetch_df(conn, COMMERCIAL_JOIN_SQL)
|
||
|
||
if joined.empty:
|
||
sections.append("\n_No contacts ↔ HRRR matches were found._\n")
|
||
else:
|
||
sections.append(
|
||
f"\n## Contact ↔ HRRR per-band Pearson correlations "
|
||
f"(matched n={len(joined)})\n"
|
||
)
|
||
sections.append(md_table(correlations_per_band(joined)))
|
||
|
||
sections.append(
|
||
"\n## Per-band HPBL bin distance distribution\n\n"
|
||
"One table per band with ≥50 matched contacts. A band-specific "
|
||
"effect here means HPBL belongs in the scorer for that band.\n"
|
||
)
|
||
for band in bands_with_samples(joined):
|
||
table = hpbl_bins(joined, band)
|
||
if table.empty:
|
||
continue
|
||
sections.append(f"\n**{band_label(band)}:**\n")
|
||
sections.append(md_table(table))
|
||
|
||
sections.append(
|
||
"\n## Per-band pressure bin distance distribution\n\n"
|
||
"Surface pressure tends to proxy synoptic-scale stagnation, so we "
|
||
"want to see longer distances under the 1015-1025 mb ridge at every "
|
||
"band that ducts.\n"
|
||
)
|
||
for band in bands_with_samples(joined):
|
||
table = pressure_bins(joined, band)
|
||
if table.empty:
|
||
continue
|
||
sections.append(f"\n**{band_label(band)}:**\n")
|
||
sections.append(md_table(table))
|
||
|
||
if joined_narr.empty:
|
||
sections.append(
|
||
"\n## Contact ↔ NARR per-band Pearson correlations\n\n"
|
||
"_No pre-2014 contacts ↔ NARR matches were found. Either the NARR "
|
||
"backfill has not run for this corpus yet or there are no "
|
||
"pre-2014 contacts with `pos1` set._\n"
|
||
)
|
||
else:
|
||
sections.append(
|
||
f"\n## Contact ↔ NARR per-band Pearson correlations "
|
||
f"(pre-2014 only, matched n={len(joined_narr)})\n\n"
|
||
"Historical sanity check: if the HRRR-era signs and magnitudes "
|
||
"survive on 30+ years of NARR reanalysis, the scoring factors are "
|
||
"physical rather than HRRR-artefact. Diverging signs = flag it.\n"
|
||
)
|
||
sections.append(md_table(correlations_per_band(joined_narr)))
|
||
|
||
sections.append(
|
||
f"\n## NEXRAD composite reflectivity vs distance "
|
||
f"(matched n={len(nexrad)})\n"
|
||
)
|
||
if nexrad.empty:
|
||
sections.append("_No contacts ↔ NEXRAD matches were found._\n")
|
||
else:
|
||
sections.append(
|
||
"Rain-attenuation sanity check — one table per band with ≥50 "
|
||
"matched contacts. At rain-sensitive bands (24+ GHz) higher "
|
||
"`max_dbz` should map to shorter contacts; at 10 GHz and below "
|
||
"the effect should be barely detectable.\n"
|
||
)
|
||
for band in bands_with_samples(nexrad):
|
||
table = nexrad_bins(nexrad, band)
|
||
if table.empty:
|
||
continue
|
||
sections.append(f"\n**{band_label(band)}:**\n")
|
||
sections.append(md_table(table))
|
||
|
||
sections.append(
|
||
f"\n## hrrr_native_profiles.best_duct_band_ghz vs distance "
|
||
f"(matched n={len(native)})\n"
|
||
)
|
||
if native.empty:
|
||
sections.append("_No contacts ↔ native-profile matches were found._\n")
|
||
else:
|
||
sections.append(
|
||
"Validates the 1.15× boost in `Scorer.score_refractivity/4`: "
|
||
"contacts where the native duct supports the target band should "
|
||
"run longer than those where it does not. One table per band with "
|
||
"≥50 matched contacts.\n"
|
||
)
|
||
for band in bands_with_samples(native):
|
||
table = native_duct_bins(native, band)
|
||
if table.empty:
|
||
continue
|
||
sections.append(f"\n**{band_label(band)}:**\n")
|
||
sections.append(md_table(table))
|
||
|
||
sections.append(
|
||
f"\n## Commercial-link rx_power degradation vs contemporaneous DFW contacts "
|
||
f"(matched n={len(commercial)})\n"
|
||
)
|
||
if commercial.empty:
|
||
sections.append(
|
||
"_No DFW-zone contacts fall inside the commercial_samples date "
|
||
"window yet. Expected to stay empty until the Aug/Sep contest "
|
||
"season produces contacts that overlap with live SNMP polling._\n"
|
||
)
|
||
else:
|
||
rho = commercial_correlation(commercial)
|
||
if rho is not None:
|
||
sections.append(f"Pearson(degradation_db, distance_km) = **{rho}**\n\n")
|
||
sections.append(md_table(commercial_degradation_summary(commercial)))
|
||
|
||
out_path.write_text("\n".join(sections))
|
||
print(f"wrote {out_path}", file=sys.stderr)
|
||
return 0
|
||
|
||
|
||
if __name__ == "__main__":
|
||
raise SystemExit(main())
|