Extend recalibrate_algo.py with NEXRAD, native-duct, commercial sections
Three new analysis passes in scripts/recalibrate_algo.py that match the
scoring factors landed in commit 86364f4:
* NEXRAD composite reflectivity vs distance per band — checks that the
Marshall-Palmer rain rate in Scorer.dbz_to_rain_rate_mmhr/1 moves the
needle in the expected direction at 24+ GHz.
* hrrr_native_profiles.best_duct_band_ghz vs distance per band —
validates the 1.15× boost in Scorer.score_refractivity/4. Current
corpus shows essentially zero contacts with a native duct supporting
≥5 GHz at the target band, consistent with the 0.12% rate in Part 2c;
the boost is physics-correct but almost never fires today.
* Commercial-link rx_power degradation vs contemporaneous DFW contacts —
pearson + binned distance. Empty today because commercial_samples only
covers Mar 30–Apr 13 and contests don't overlap. Scaffolding is in
place so the moment contest-season overlap exists, the signal shows.
Report rebuilt at docs/algo-reports/2026-04-13-recalibration.md.
This commit is contained in:
parent
fa26db897d
commit
5af4142e28
2 changed files with 324 additions and 13 deletions
|
|
@ -15,13 +15,13 @@ Statement timeout: 20min
|
|||
| contacts | 58560 | 1991-05-04 | 2026-03-07 |
|
||||
| era5_profiles | 0 | | |
|
||||
| hrrr_native_profiles | 11472 | 2019-03-09 | 2025-09-21 |
|
||||
| hrrr_profiles | 76058705 | 2016-09-17 | 2026-04-13 |
|
||||
| hrrr_profiles | 76151193 | 2016-09-17 | 2026-04-13 |
|
||||
| iemre_observations | 15745 | | |
|
||||
| nexrad_observations | 7286 | 2019-08-17 | 2024-09-22 |
|
||||
| propagation_scores | 9618752 | 2026-04-13 | 2026-04-13 |
|
||||
| propagation_scores | 10543632 | 2026-04-13 | 2026-04-13 |
|
||||
| rtma_observations | 0 | | |
|
||||
| soundings | 16864 | 1990-06-07 | 2026-03-07 |
|
||||
| surface_observations | 392442 | 1991-05-03 | 2026-04-13 |
|
||||
| surface_observations | 392444 | 1991-05-03 | 2026-04-13 |
|
||||
| terrain_profiles | 58554 | | |
|
||||
|
||||
|
||||
|
|
@ -93,20 +93,20 @@ Statement timeout: 20min
|
|||
| band_mhz | n | rho_pr | rho_dpc | rho_pwat | rho_sref | rho_grad | rho_tc | rho_hpbl |
|
||||
|-----------:|--------:|---------:|----------:|-----------:|-----------:|-----------:|---------:|-----------:|
|
||||
| 1296.000 | 63.000 | -0.446 | 0.150 | -0.004 | 0.321 | -0.308 | -0.313 | -0.383 |
|
||||
| 10000.000 | 680.000 | -0.104 | 0.046 | 0.003 | 0.181 | -0.127 | -0.108 | -0.200 |
|
||||
| 24000.000 | 147.000 | -0.290 | -0.179 | -0.220 | 0.096 | -0.034 | -0.301 | -0.226 |
|
||||
| 47000.000 | 75.000 | -0.111 | 0.015 | 0.006 | 0.218 | -0.167 | 0.086 | -0.154 |
|
||||
| 10000.000 | 680.000 | -0.109 | 0.044 | 0.002 | 0.182 | -0.127 | -0.107 | -0.198 |
|
||||
| 24000.000 | 147.000 | -0.290 | -0.179 | -0.219 | 0.096 | -0.036 | -0.300 | -0.223 |
|
||||
| 47000.000 | 75.000 | -0.111 | 0.016 | 0.007 | 0.217 | -0.175 | 0.089 | -0.154 |
|
||||
|
||||
|
||||
## 10 GHz: HPBL bin distance distribution
|
||||
|
||||
| bin | n | avg_km | p50_km |
|
||||
|:----------|----:|---------:|---------:|
|
||||
| <200 | 55 | 222.600 | 174.000 |
|
||||
| 200-500 | 174 | 183.200 | 154.500 |
|
||||
| 500-1000 | 181 | 163.700 | 126.000 |
|
||||
| 1000-1500 | 103 | 160.900 | 136.000 |
|
||||
| 1500-2000 | 93 | 140.500 | 128.000 |
|
||||
| <200 | 54 | 217.600 | 174.000 |
|
||||
| 200-500 | 177 | 184.800 | 160.000 |
|
||||
| 500-1000 | 177 | 163.800 | 126.000 |
|
||||
| 1000-1500 | 104 | 158.800 | 136.000 |
|
||||
| 1500-2000 | 94 | 142.800 | 136.000 |
|
||||
| >=2000 | 74 | 100.300 | 90.000 |
|
||||
|
||||
|
||||
|
|
@ -114,8 +114,67 @@ Statement timeout: 20min
|
|||
|
||||
| bin | n | avg_km | p50_km |
|
||||
|:----------|----:|---------:|---------:|
|
||||
| <990 | 101 | 211.900 | 134.000 |
|
||||
| 990-1000 | 538 | 155.400 | 128.000 |
|
||||
| <990 | 103 | 214.200 | 136.000 |
|
||||
| 990-1000 | 536 | 154.800 | 128.000 |
|
||||
| 1000-1010 | 28 | 164.600 | 115.000 |
|
||||
| 1010-1020 | 2 | 23.000 | 23.000 |
|
||||
| >=1020 | 11 | 103.100 | 23.000 |
|
||||
|
||||
|
||||
## NEXRAD composite reflectivity vs distance (matched n=102)
|
||||
|
||||
Sanity check for the rain-attenuation direction: at rain-sensitive bands higher `max_dbz` should map to shorter contacts. At 10 GHz the effect should be barely detectable.
|
||||
|
||||
|
||||
**10 GHz:**
|
||||
|
||||
| bin | n | avg_km | p50_km |
|
||||
|:-----------------|----:|---------:|---------:|
|
||||
| none (<5) | 59 | 155.700 | 143.000 |
|
||||
| drizzle (5-20) | 7 | 56.400 | 49.000 |
|
||||
| light (20-30) | 13 | 95.200 | 40.000 |
|
||||
| moderate (30-45) | 4 | 97.200 | 92.000 |
|
||||
| heavy (>=45) | 1 | 213.000 | 213.000 |
|
||||
|
||||
|
||||
**24 GHz:**
|
||||
|
||||
| bin | n | avg_km | p50_km |
|
||||
|:-----------------|----:|---------:|---------:|
|
||||
| none (<5) | 5 | 29.800 | 23.000 |
|
||||
| drizzle (5-20) | 1 | 49.000 | 49.000 |
|
||||
| light (20-30) | 5 | 34.000 | 40.000 |
|
||||
| moderate (30-45) | 1 | 60.000 | 60.000 |
|
||||
|
||||
|
||||
## hrrr_native_profiles.best_duct_band_ghz vs distance (matched n=619)
|
||||
|
||||
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.
|
||||
|
||||
|
||||
**10 GHz:**
|
||||
|
||||
| bin | n | avg_km | p50_km |
|
||||
|:-------|--------:|---------:|---------:|
|
||||
| none | 447.000 | 140.800 | 128.000 |
|
||||
| <5 GHz | 38.000 | 150.000 | 114.500 |
|
||||
|
||||
|
||||
**24 GHz:**
|
||||
|
||||
| bin | n | avg_km | p50_km |
|
||||
|:-------|-------:|---------:|---------:|
|
||||
| none | 97.000 | 47.600 | 32.000 |
|
||||
| <5 GHz | 3.000 | 56.300 | 60.000 |
|
||||
|
||||
|
||||
**47 GHz:**
|
||||
|
||||
| bin | n | avg_km | p50_km |
|
||||
|:------|-------:|---------:|---------:|
|
||||
| none | 16.000 | 54.700 | 42.000 |
|
||||
|
||||
|
||||
## Commercial-link rx_power degradation vs contemporaneous DFW contacts (matched n=0)
|
||||
|
||||
_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._
|
||||
|
|
|
|||
|
|
@ -12,6 +12,9 @@ and produces a Markdown report covering:
|
|||
* native HRRR profile duct stats
|
||||
* per-band Pearson correlations of contact distance vs HRRR fields
|
||||
* binned distance distributions for HPBL and pressure at 10 GHz
|
||||
* 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.
|
||||
|
|
@ -148,6 +151,102 @@ WITH joined AS (
|
|||
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.
|
||||
DATA_GAP_SQL = """
|
||||
SELECT
|
||||
|
|
@ -226,6 +325,98 @@ def pressure_bins(df: pd.DataFrame, band_mhz: int = 10000) -> pd.DataFrame:
|
|||
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 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
|
||||
|
|
@ -320,6 +511,15 @@ def main() -> int:
|
|||
print("• per-band correlations (this is the slow one)", file=sys.stderr)
|
||||
joined = fetch_df(conn, PER_BAND_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:
|
||||
|
|
@ -335,6 +535,58 @@ def main() -> int:
|
|||
sections.append("\n## 10 GHz: pressure bin distance distribution\n")
|
||||
sections.append(md_table(pressure_bins(joined)))
|
||||
|
||||
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(
|
||||
"Sanity check for the rain-attenuation direction: at rain-sensitive "
|
||||
"bands higher `max_dbz` should map to shorter contacts. At 10 GHz "
|
||||
"the effect should be barely detectable.\n\n"
|
||||
)
|
||||
sections.append("**10 GHz:**\n")
|
||||
sections.append(md_table(nexrad_bins(nexrad, 10_000)))
|
||||
sections.append("\n**24 GHz:**\n")
|
||||
sections.append(md_table(nexrad_bins(nexrad, 24_000)))
|
||||
|
||||
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.\n\n"
|
||||
)
|
||||
sections.append("**10 GHz:**\n")
|
||||
sections.append(md_table(native_duct_bins(native, 10_000)))
|
||||
sections.append("\n**24 GHz:**\n")
|
||||
sections.append(md_table(native_duct_bins(native, 24_000)))
|
||||
sections.append("\n**47 GHz:**\n")
|
||||
sections.append(md_table(native_duct_bins(native, 47_000)))
|
||||
|
||||
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
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue