Graham McIntire
1086f52c85
Both RTMA and METAR5 schemas + clients + workers were defined but never had a consumer. `find_nearest_rtma/3` and `recent_surface_obs/3` existed in `Microwaveprop.Weather` with zero callers cluster-wide; RtmaFetchWorker had test coverage but was never enqueued anywhere; no IEM 5-min ASOS fetcher was ever written. The tables sat empty in prod and the recalibrate audit was permanently flagging them BROKEN even though the empty state was the steady state. Drop the lot — schemas, clients, workers, queue config, Req.Test stubs, audit checks, and the per-table unique tests in the observation-changesets suite. New migration drops the now-unreferenced `rtma_observations` and `metar_5min_observations` tables (both 0 rows in prod). Trim the `:rtma` slot out of `backfill_only_queues` in runtime.exs so Oban Pro's Smart engine stops trying to manage a queue with no producer. Audit thresholds reset on the surviving NARR check: drop the "WARN < 5000 absolute rows" rule that didn't account for NarrFetchWorker's 0.13° spatial dedup, and point remediation at the existing BackfillEnqueueWorker cron instead of the dev-only `mix narr.backfill` task. scripts/recalibrate_algo.py: ROW_COUNT_SOURCES + DATA_GAP_SQL + DATA_GAP_RULES all stop referencing the dropped tables.
916 lines
35 KiB
Python
Executable file
916 lines
35 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"),
|
||
("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 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",
|
||
# NARR feeds pre-2014 enrichment via BackfillEnqueueWorker's
|
||
# 30-min cron. NarrFetchWorker dedupes by 0.13° grid, so the
|
||
# equilibrium count is ~3 path points × pre2014_contacts after
|
||
# spatial collapse — not the contact count itself. Threshold
|
||
# follows: BROKEN at zero, WARN below half the expected
|
||
# post-dedup ceiling.
|
||
"broken_when": lambda v: v == 0,
|
||
"warn_when": lambda _v: False,
|
||
"remediation": (
|
||
"BackfillEnqueueWorker dispatches narr jobs every 30 min; "
|
||
"if rows aren't climbing, check NarrFetchWorker for upstream "
|
||
"NCEI errors in the Oban dashboard."
|
||
),
|
||
},
|
||
{
|
||
"key": "narr_pre2014_rows",
|
||
"broken_when": lambda v: v == 0,
|
||
"warn_when": lambda _v: False,
|
||
"remediation": (
|
||
"NARR's purpose is filling the pre-2014-10 gap; without rows here "
|
||
"the historical-sanity correlation table is empty. Same self-heal "
|
||
"as narr_rows."
|
||
),
|
||
},
|
||
{
|
||
"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": "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())
|