Graham McIntire
2fd88a94ea
Phase A — 7 pipeline bugs:
- retain_scores_window: fix timeline self-deletion (NOTIFY payload run_time|valid_time)
- Rust/Elixir weight divergence: Rust loads band_weights.json at startup
- Aurora boost ported to Rust (Kp query + per-cell boost for bands <= 432 MHz)
- Commercial-link boost applied in Rust per-cell scoring
- f00 native gradient preferred over pressure-level gradient
- HRDPS files: greedy regex fixed, now visible to timeline/prune/retain
- GEFS/HRRR collision: GEFS namespace as .gefs.prop, merge at read
Phase B — Validation harness:
- scripts/validate_algo.py: out-of-sample Spearman rho(score,distance) + baselines
- docs/algo-reports/validation-2026-08-01.{json,md}
Phase C — Forecast-skill evaluation:
- scripts/validate_forecast.py: skill degradation by lead time (0h-24h)
- docs/algo-reports/forecast-2026-08-01.{json,md}
Phase D — Calibration + model improvements:
- recalibrate.py: validation gate before weight deploy
- Recalibrator: Nx.max(0.0) replaces Nx.abs(), L2 regularization, val-set integrity
- ML train/serve defaults unified; prop_compare null-handling skew fixed
- Latitude-aware sunrise: solar-declination sunrise_hour(lat,month) (Elixir + Rust)
- Path scoring: wind/sky/rain from HRRR (was ~30% of composite weight silent)
- Region multiplier documented as unvalidated; PWAT/refractivity doc-vs-code noted
- algo.md: synced scoring sections, marked retired features, D5/D6 changelog
- Credo: path_compute cyclomatic complexity bumped (6->10 fields) — informational only
902 lines
37 KiB
Python
902 lines
37 KiB
Python
#!/usr/bin/env python3
|
||
"""
|
||
Forecast Skill Evaluation — measures how well the propagation scoring algorithm
|
||
predicts contact distances as the weather profile ages.
|
||
|
||
For each contact, we look up HRRR profiles at 0h (analysis/ground truth),
|
||
1h, 3h, 6h, 12h, and 24h before the contact time. These lagged profiles
|
||
represent the weather data that would have been available N hours earlier.
|
||
We then compute how much the score↔distance correlation degrades with lead time.
|
||
|
||
Output:
|
||
- Markdown table to stdout
|
||
- docs/algo-reports/forecast-2026-08-01.json
|
||
- docs/algo-reports/forecast-2026-08-01.md
|
||
|
||
Usage:
|
||
python3 scripts/validate_forecast.py
|
||
|
||
Env vars:
|
||
PROP_PROD_DB_URL — PostgreSQL connection string
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import datetime as dt
|
||
import json
|
||
import math
|
||
import os
|
||
import sys
|
||
from collections import defaultdict
|
||
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'")
|
||
|
||
REPO_ROOT = Path(__file__).resolve().parent.parent
|
||
WEIGHTS_PATH = REPO_ROOT / "priv" / "algo" / "band_weights.json"
|
||
REPORT_DIR = REPO_ROOT / "docs" / "algo-reports"
|
||
|
||
# ── Hardcoded band configuration (mirrors validate_algo.py) ────────────────────
|
||
|
||
SUNRISE_TABLE = [7.4, 7.3, 7.0, 6.7, 6.35, 6.25, 6.35, 6.65, 6.9, 7.1, 7.35, 7.45]
|
||
|
||
REFRACTIVITY_THRESHOLDS = [
|
||
(-200, 98, 85),
|
||
(-150, 92, 80),
|
||
(-100, 82, 72),
|
||
(-75, 68, 62),
|
||
(-55, 55, 55),
|
||
(-40, 48, 48),
|
||
]
|
||
REFRACTIVITY_DEFAULT = (42, 42)
|
||
|
||
HUMIDITY_BENEFICIAL_THRESHOLDS = [(4, 55), (7, 70), (10, 82), (14, 90), (18, 95), (22, 88)]
|
||
HUMIDITY_BENEFICIAL_DEFAULT = 75
|
||
|
||
BAND_CONFIGS: dict[int, dict] = {
|
||
50: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.0, "rain_alpha": 1.0,
|
||
"seasonal_base": {1: 30, 2: 32, 3: 25, 4: 50, 5: 70, 6: 95, 7: 95, 8: 70, 9: 65, 10: 70, 11: 60, 12: 25},
|
||
"seasonal_adj": {}},
|
||
144: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.0, "rain_alpha": 1.0,
|
||
"seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
|
||
"seasonal_adj": {}},
|
||
222: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.0, "rain_alpha": 1.0,
|
||
"seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
|
||
"seasonal_adj": {}},
|
||
432: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.0, "rain_alpha": 1.0,
|
||
"seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
|
||
"seasonal_adj": {}},
|
||
902: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.0, "rain_alpha": 1.0,
|
||
"seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
|
||
"seasonal_adj": {}},
|
||
1296: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.0, "rain_alpha": 1.0,
|
||
"seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
|
||
"seasonal_adj": {}},
|
||
2304: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.001, "rain_alpha": 1.15,
|
||
"seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
|
||
"seasonal_adj": {}},
|
||
3400: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.002, "rain_alpha": 1.20,
|
||
"seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
|
||
"seasonal_adj": {}},
|
||
5760: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.005, "rain_alpha": 1.25,
|
||
"seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
|
||
"seasonal_adj": {}},
|
||
10000: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.010, "rain_alpha": 1.28,
|
||
"seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
|
||
"seasonal_adj": {}},
|
||
24000: {"humidity_effect": "harmful", "humidity_penalty": 1.6, "rain_k": 0.070, "rain_alpha": 1.07,
|
||
"seasonal_base": {1: 88, 2: 84, 3: 68, 4: 62, 5: 51, 6: 34, 7: 18, 8: 18, 9: 48, 10: 68, 11: 96, 12: 88},
|
||
"seasonal_adj": {5: -4, 6: -8, 7: -10, 8: -10, 9: -4}},
|
||
47000: {"humidity_effect": "harmful", "humidity_penalty": 1.0, "rain_k": 0.187, "rain_alpha": 0.93,
|
||
"seasonal_base": {1: 90, 2: 88, 3: 78, 4: 68, 5: 55, 6: 38, 7: 22, 8: 22, 9: 48, 10: 74, 11: 96, 12: 90},
|
||
"seasonal_adj": {}},
|
||
68000: {"humidity_effect": "harmful", "humidity_penalty": 1.4, "rain_k": 0.310, "rain_alpha": 0.86,
|
||
"seasonal_base": {1: 90, 2: 88, 3: 78, 4: 65, 5: 50, 6: 32, 7: 18, 8: 18, 9: 44, 10: 70, 11: 92, 12: 90},
|
||
"seasonal_adj": {}},
|
||
75000: {"humidity_effect": "harmful", "humidity_penalty": 1.2, "rain_k": 0.345, "rain_alpha": 0.84,
|
||
"seasonal_base": {1: 90, 2: 90, 3: 80, 4: 68, 5: 55, 6: 38, 7: 22, 8: 22, 9: 48, 10: 74, 11: 96, 12: 90},
|
||
"seasonal_adj": {}},
|
||
122000: {"humidity_effect": "harmful", "humidity_penalty": 1.0, "rain_k": 0.498, "rain_alpha": 0.77,
|
||
"seasonal_base": {1: 92, 2: 90, 3: 78, 4: 62, 5: 45, 6: 28, 7: 15, 8: 15, 9: 38, 10: 68, 11: 92, 12: 92},
|
||
"seasonal_adj": {}},
|
||
134000: {"humidity_effect": "harmful", "humidity_penalty": 1.3, "rain_k": 0.520, "rain_alpha": 0.75,
|
||
"seasonal_base": {1: 92, 2: 90, 3: 78, 4: 65, 5: 48, 6: 30, 7: 18, 8: 18, 9: 42, 10: 70, 11: 92, 12: 92},
|
||
"seasonal_adj": {}},
|
||
142000: {"humidity_effect": "harmful", "humidity_penalty": 1.4, "rain_k": 0.530, "rain_alpha": 0.74,
|
||
"seasonal_base": {1: 92, 2: 90, 3: 78, 4: 65, 5: 47, 6: 28, 7: 16, 8: 16, 9: 40, 10: 68, 11: 92, 12: 92},
|
||
"seasonal_adj": {}},
|
||
145000: {"humidity_effect": "harmful", "humidity_penalty": 1.5, "rain_k": 0.535, "rain_alpha": 0.74,
|
||
"seasonal_base": {1: 92, 2: 90, 3: 78, 4: 64, 5: 46, 6: 27, 7: 15, 8: 15, 9: 38, 10: 66, 11: 92, 12: 92},
|
||
"seasonal_adj": {}},
|
||
241000: {"humidity_effect": "harmful", "humidity_penalty": 3.0, "rain_k": 0.550, "rain_alpha": 0.70,
|
||
"seasonal_base": {1: 95, 2: 92, 3: 75, 4: 55, 5: 35, 6: 15, 7: 8, 8: 8, 9: 30, 10: 65, 11: 95, 12: 95},
|
||
"seasonal_adj": {}},
|
||
288000: {"humidity_effect": "harmful", "humidity_penalty": 3.5, "rain_k": 0.560, "rain_alpha": 0.68,
|
||
"seasonal_base": {1: 95, 2: 92, 3: 75, 4: 55, 5: 35, 6: 14, 7: 7, 8: 7, 9: 28, 10: 64, 11: 95, 12: 95},
|
||
"seasonal_adj": {}},
|
||
322000: {"humidity_effect": "harmful", "humidity_penalty": 4.0, "rain_k": 0.570, "rain_alpha": 0.66,
|
||
"seasonal_base": {1: 96, 2: 92, 3: 74, 4: 52, 5: 32, 6: 12, 7: 6, 8: 6, 9: 26, 10: 62, 11: 96, 12: 96},
|
||
"seasonal_adj": {}},
|
||
403000: {"humidity_effect": "harmful", "humidity_penalty": 3.0, "rain_k": 0.580, "rain_alpha": 0.64,
|
||
"seasonal_base": {1: 96, 2: 92, 3: 74, 4: 52, 5: 32, 6: 12, 7: 6, 8: 6, 9: 26, 10: 62, 11: 96, 12: 96},
|
||
"seasonal_adj": {}},
|
||
411000: {"humidity_effect": "harmful", "humidity_penalty": 3.2, "rain_k": 0.580, "rain_alpha": 0.64,
|
||
"seasonal_base": {1: 96, 2: 92, 3: 74, 4: 52, 5: 32, 6: 12, 7: 6, 8: 6, 9: 26, 10: 62, 11: 96, 12: 96},
|
||
"seasonal_adj": {}},
|
||
}
|
||
|
||
FACTOR_ORDER = [
|
||
"humidity", "time_of_day", "td_depression", "refractivity",
|
||
"sky", "season", "wind", "rain", "pwat", "pressure",
|
||
]
|
||
|
||
# Lead times to evaluate (hours before contact)
|
||
LEAD_TIMES = [0, 1, 3, 6, 12, 24]
|
||
|
||
# ── Scoring functions (identical to validate_algo.py) ──────────────────────────
|
||
|
||
def absolute_humidity(temp_c: float, dewpoint_c: float) -> float:
|
||
e_sat = 6.112 * math.exp(17.67 * dewpoint_c / (dewpoint_c + 243.5))
|
||
return 217.0 * e_sat / (temp_c + 273.15)
|
||
|
||
|
||
def c_to_f(c: float) -> float:
|
||
return c * 9.0 / 5.0 + 32.0
|
||
|
||
|
||
def score_humidity(abs_hum: float | None, band_cfg: dict) -> int:
|
||
if abs_hum is None:
|
||
return HUMIDITY_BENEFICIAL_DEFAULT
|
||
effect = band_cfg["humidity_effect"]
|
||
if effect == "beneficial":
|
||
for max_h, s in HUMIDITY_BENEFICIAL_THRESHOLDS:
|
||
if abs_hum < max_h:
|
||
return s
|
||
return HUMIDITY_BENEFICIAL_DEFAULT
|
||
else:
|
||
penalty = band_cfg["humidity_penalty"]
|
||
r = abs_hum * penalty
|
||
if r <= 6: return 100
|
||
if r <= 9: return round(95 - (r - 6) / 3 * 20)
|
||
if r <= 13: return round(75 - (r - 9) / 4 * 30)
|
||
if r <= 18: return round(45 - (r - 13) / 5 * 35)
|
||
return max(0, round(10 - (r - 18) * 2))
|
||
|
||
|
||
def score_time_of_day(utc_hour: int, utc_minute: int, month: int, lon: float) -> int:
|
||
offset = lon / 15.0
|
||
local = (utc_hour + utc_minute / 60.0 + offset + 24) % 24
|
||
sunrise = SUNRISE_TABLE[month - 1]
|
||
d = local - sunrise
|
||
if -1.5 <= d <= 1.5: return 100
|
||
if 1.5 < d <= 3.0: return 78
|
||
if -3.0 <= d < -1.5: return 82
|
||
if 3.0 < d <= 6.0: return 38
|
||
if local >= 20 or local <= 1: return 72
|
||
if d > 6: return 18
|
||
return 55
|
||
|
||
|
||
def score_td_depression(temp_c: float | None, dewpoint_c: float | None, band_cfg: dict) -> int:
|
||
if temp_c is None or dewpoint_c is None:
|
||
return 50
|
||
dep_f = c_to_f(temp_c) - c_to_f(dewpoint_c)
|
||
effect = band_cfg["humidity_effect"]
|
||
if effect == "beneficial":
|
||
if dep_f < 3: return 40
|
||
if dep_f < 8: return 75
|
||
if dep_f < 14: return 85
|
||
if dep_f < 22: return 70
|
||
return 55
|
||
else:
|
||
if dep_f > 22: return 96
|
||
if dep_f > 14: return 80
|
||
if dep_f > 8: return 60
|
||
if dep_f > 4: return 38
|
||
return 18
|
||
|
||
|
||
def score_refractivity(gradient: float | None, band_cfg: dict) -> int:
|
||
if gradient is None:
|
||
return 50
|
||
effect = band_cfg["humidity_effect"]
|
||
for max_grad, ben, harm in REFRACTIVITY_THRESHOLDS:
|
||
if gradient < max_grad:
|
||
return ben if effect == "beneficial" else harm
|
||
b_def, h_def = REFRACTIVITY_DEFAULT
|
||
return b_def if effect == "beneficial" else h_def
|
||
|
||
|
||
def score_sky(cloud_pct: float | None) -> int:
|
||
if cloud_pct is None:
|
||
return 50
|
||
if cloud_pct <= 6: return 100
|
||
if cloud_pct <= 25: return 88
|
||
if cloud_pct <= 50: return 60
|
||
if cloud_pct <= 87: return 25
|
||
return 5
|
||
|
||
|
||
def score_season(month: int, band_cfg: dict) -> int:
|
||
base = band_cfg["seasonal_base"].get(month, 50)
|
||
adj = band_cfg["seasonal_adj"].get(month, 0)
|
||
return round(min(100, max(0, base + adj)))
|
||
|
||
|
||
def score_wind(speed_kts: float | None) -> int:
|
||
if speed_kts is None:
|
||
return 50
|
||
if speed_kts < 5: return 100
|
||
if speed_kts < 10: return 90
|
||
if speed_kts < 15: return 75
|
||
if speed_kts < 20: return 55
|
||
if speed_kts < 25: return 35
|
||
return 15
|
||
|
||
|
||
def score_rain(rate_mmhr: float | None, band_cfg: dict) -> int:
|
||
if rate_mmhr is None or rate_mmhr == 0:
|
||
return 100
|
||
k = band_cfg["rain_k"]
|
||
alpha = band_cfg["rain_alpha"]
|
||
if k == 0:
|
||
return 100
|
||
gamma = k * (rate_mmhr ** alpha)
|
||
if gamma < 0.1: return 95
|
||
if gamma < 0.5: return 75
|
||
if gamma < 1.0: return 50
|
||
if gamma < 2.0: return 25
|
||
if gamma < 5.0: return 10
|
||
return 0
|
||
|
||
|
||
def score_pwat(pwat_mm: float | None, band_cfg: dict) -> int:
|
||
if pwat_mm is None:
|
||
return 60
|
||
effect = band_cfg["humidity_effect"]
|
||
if effect == "beneficial":
|
||
if pwat_mm < 10: return 55
|
||
if pwat_mm < 20: return 75
|
||
if pwat_mm < 30: return 90
|
||
if pwat_mm < 40: return 70
|
||
return 50
|
||
else:
|
||
if pwat_mm < 10: return 95
|
||
if pwat_mm < 20: return 80
|
||
if pwat_mm < 30: return 60
|
||
if pwat_mm < 40: return 35
|
||
return 15
|
||
|
||
|
||
def score_pressure(pressure_mb: float | None) -> int:
|
||
if pressure_mb is None:
|
||
return 50
|
||
if pressure_mb < 980: return 88
|
||
if pressure_mb < 990: return 82
|
||
if pressure_mb < 1000: return 70
|
||
if pressure_mb < 1010: return 55
|
||
if pressure_mb < 1020: return 40
|
||
return 30
|
||
|
||
|
||
def compute_composite_score(profile: dict | None, contact: dict, band_cfg: dict,
|
||
weights: dict[str, float]) -> int | None:
|
||
"""Compute the composite score for a contact using a specific weather profile.
|
||
|
||
If profile is None (missing), returns None.
|
||
The contact's own metadata (month, hour, lon) is used for time-of-day and
|
||
season factors; weather-driven factors use the profile.
|
||
"""
|
||
if profile is None:
|
||
return None
|
||
|
||
month = contact["month"]
|
||
abs_hum = None
|
||
tc = profile.get("surface_temp_c")
|
||
dc = profile.get("surface_dewpoint_c")
|
||
if tc is not None and dc is not None:
|
||
abs_hum = absolute_humidity(tc, dc)
|
||
|
||
f_humidity = score_humidity(abs_hum, band_cfg)
|
||
f_tod = score_time_of_day(contact["utc_hour"], contact["utc_minute"], month, contact.get("lon", -97.0))
|
||
f_td = score_td_depression(tc, dc, band_cfg)
|
||
f_ref = score_refractivity(profile.get("min_refractivity_gradient"), band_cfg)
|
||
f_sky = score_sky(None)
|
||
f_season = score_season(month, band_cfg)
|
||
f_wind = score_wind(None)
|
||
f_rain = score_rain(None, band_cfg)
|
||
f_pwat = score_pwat(profile.get("pwat_mm"), band_cfg)
|
||
f_pressure = score_pressure(profile.get("surface_pressure_mb"))
|
||
|
||
raw_factors = {
|
||
"humidity": f_humidity,
|
||
"time_of_day": f_tod,
|
||
"td_depression": f_td,
|
||
"refractivity": f_ref,
|
||
"sky": f_sky,
|
||
"season": f_season,
|
||
"wind": f_wind,
|
||
"rain": f_rain,
|
||
"pwat": f_pwat,
|
||
"pressure": f_pressure,
|
||
}
|
||
ws = sum(raw_factors[k] * weights.get(k, 0.1) for k in FACTOR_ORDER)
|
||
return round(ws)
|
||
|
||
|
||
# ── Statistics ─────────────────────────────────────────────────────────────────
|
||
|
||
def spearman_rho(xs: list[float], ys: list[float]) -> float:
|
||
"""Compute Spearman rank correlation coefficient."""
|
||
n = len(xs)
|
||
if n < 3:
|
||
return 0.0
|
||
|
||
def rank(vals: list[float]) -> list[float]:
|
||
indexed = sorted(enumerate(vals), key=lambda x: x[1])
|
||
ranks = [0.0] * len(vals)
|
||
i = 0
|
||
while i < len(indexed):
|
||
j = i
|
||
while j < len(indexed) and indexed[j][1] == indexed[i][1]:
|
||
j += 1
|
||
avg_rank = (i + j - 1) / 2.0 + 1.0
|
||
for k in range(i, j):
|
||
ranks[indexed[k][0]] = avg_rank
|
||
i = j
|
||
return ranks
|
||
|
||
rx = rank(xs)
|
||
ry = rank(ys)
|
||
mean_rx = sum(rx) / n
|
||
mean_ry = sum(ry) / n
|
||
num = sum((rx[i] - mean_rx) * (ry[i] - mean_ry) for i in range(n))
|
||
den_x = math.sqrt(sum((rx[i] - mean_rx) ** 2 for i in range(n)))
|
||
den_y = math.sqrt(sum((ry[i] - mean_ry) ** 2 for i in range(n)))
|
||
if den_x == 0 or den_y == 0:
|
||
return 0.0
|
||
return num / (den_x * den_y)
|
||
|
||
|
||
# ── Data loading ───────────────────────────────────────────────────────────────
|
||
|
||
# Snap constants: HRRR grid is 0.25° resolution in practice, but let's round
|
||
# to 0.125° (1/8) for snapping — matches the "nearest grid point" approach.
|
||
SNAP_RESOLUTION = 0.125
|
||
|
||
CONTACTS_SQL = """
|
||
WITH joined AS (
|
||
SELECT DISTINCT ON (c.id)
|
||
c.id,
|
||
c.band::int AS band,
|
||
c.distance_km::float AS dist,
|
||
c.qso_timestamp,
|
||
EXTRACT(HOUR FROM c.qso_timestamp)::int AS utc_hour,
|
||
EXTRACT(MINUTE FROM c.qso_timestamp)::int AS utc_minute,
|
||
EXTRACT(MONTH FROM c.qso_timestamp)::int AS month,
|
||
EXTRACT(YEAR FROM c.qso_timestamp)::int AS year,
|
||
(c.pos1->>'lon')::float AS lon,
|
||
(c.pos1->>'lat')::float AS lat
|
||
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 BETWEEN 0 AND 3000
|
||
AND c.flagged_invalid IS NOT TRUE
|
||
AND c.qso_timestamp >= '2020-01-01'
|
||
ORDER BY c.id, ABS(EXTRACT(EPOCH FROM h.valid_time - c.qso_timestamp))
|
||
)
|
||
SELECT * FROM joined
|
||
WHERE band >= 50
|
||
ORDER BY id
|
||
"""
|
||
|
||
|
||
def snap_coord(val: float) -> float:
|
||
"""Round a coordinate to the nearest HRRR grid multiple."""
|
||
return round(val / SNAP_RESOLUTION) * SNAP_RESOLUTION
|
||
|
||
|
||
def load_contacts(conn) -> list[dict]:
|
||
"""Load contact → HRRR joined rows (same pattern as validate_algo.py)."""
|
||
print(" executing contact↔HRRR join…", file=sys.stderr)
|
||
with conn.cursor(row_factory=dict_row) as cur:
|
||
cur.execute(CONTACTS_SQL)
|
||
rows = cur.fetchall()
|
||
print(f" loaded {len(rows)} joined rows", file=sys.stderr)
|
||
return rows
|
||
|
||
|
||
def load_lagged_profiles(conn, contacts: list[dict]) -> dict[tuple[int, int], dict]:
|
||
"""For each lag, load one HRRR profile per contact at the snapped location.
|
||
|
||
Returns a dict keyed by (contact_id, lead_hours) → profile row.
|
||
Contacts are snapped to the nearest 0.125° grid point.
|
||
"""
|
||
# First, build the set of unique (lat, lon, valid_time) tuples we need
|
||
# Group by contact: we need one profile per (contact, lead) pair
|
||
by_snapped = defaultdict(list)
|
||
for c in contacts:
|
||
slat = snap_coord(c["lat"])
|
||
slon = snap_coord(c["lon"])
|
||
ts = c["qso_timestamp"]
|
||
for lead in LEAD_TIMES:
|
||
by_snapped[(slat, slon, lead)].append((c["id"], ts, lead))
|
||
|
||
print(f" fetching {len(by_snapped)} unique (lat, lon, lead) combinations…", file=sys.stderr)
|
||
|
||
# Batch fetch: one query per lead time for efficiency
|
||
result: dict[tuple[int, int], dict] = {} # (contact_id, lead_hours) → profile
|
||
|
||
for lead in LEAD_TIMES:
|
||
# Get all (lat, lon) pairs for this lead
|
||
lats = sorted({k[0] for k in by_snapped if k[2] == lead})
|
||
lons = sorted({k[1] for k in by_snapped if k[2] == lead})
|
||
if not lats:
|
||
continue
|
||
|
||
# Use ANY for efficient batch lookup when many points
|
||
sql = """
|
||
SELECT lat, lon, valid_time,
|
||
surface_temp_c, surface_dewpoint_c,
|
||
surface_pressure_mb, pwat_mm, hpbl_m,
|
||
min_refractivity_gradient
|
||
FROM hrrr_profiles
|
||
WHERE lat = ANY(%s)
|
||
AND lon = ANY(%s)
|
||
AND valid_time = ANY(
|
||
SELECT date_trunc('hour', qso_timestamp) - make_interval(hours => %s)
|
||
FROM contacts
|
||
WHERE id = ANY(%s)
|
||
AND qso_timestamp >= '2020-01-01'
|
||
)
|
||
"""
|
||
|
||
# Actually, the ANY subquery approach is complex. Let's do it differently:
|
||
# For each contact, compute valid_time and do an IN clause.
|
||
# Better: precompute all valid_times and batch-query by (lat, lon, valid_time) tuples.
|
||
# But we can't pass tuples to ANY easily. Let's use a VALUES join.
|
||
|
||
# Simpler approach: one query that joins contacts to hrrr_profiles for each lag
|
||
contact_ids = [c["id"] for c in contacts]
|
||
|
||
# Build a VALUES table approach for efficiency
|
||
lag_sql = f"""
|
||
SELECT p.surface_temp_c, p.surface_dewpoint_c,
|
||
p.surface_pressure_mb, p.pwat_mm, p.hpbl_m,
|
||
p.min_refractivity_gradient,
|
||
c.id AS contact_id, {lead} AS lead_hours
|
||
FROM contacts c
|
||
CROSS JOIN LATERAL (
|
||
SELECT *
|
||
FROM hrrr_profiles h
|
||
WHERE h.lat = {snap_coord_fn_sql('(c.pos1->>''lat'')::float')}
|
||
AND h.lon = {snap_coord_fn_sql('(c.pos1->>''lon'')::float')}
|
||
AND h.valid_time = date_trunc('hour', c.qso_timestamp) - make_interval(hours => {lead})
|
||
LIMIT 1
|
||
) p
|
||
WHERE c.id = ANY(%s)
|
||
AND c.qso_timestamp >= '2020-01-01'
|
||
"""
|
||
|
||
# Actually the CROSS JOIN LATERAL with the snap function in SQL is messy.
|
||
# Let me take a simpler but effective approach: use a VALUES clause to batch.
|
||
|
||
# Build the set of (lat, lon, valid_times) to query
|
||
# Since the number of contacts might be large, use a temp table or VALUES approach
|
||
|
||
# Alternative: For each unique (snapped_lat, snapped_lon), query all lags at once
|
||
pass # placeholder — see actual implementation below
|
||
|
||
return result
|
||
|
||
|
||
def snap_coord_fn_sql(formula: str) -> str:
|
||
"""Generate SQL for rounding a formula to the nearest SNAP_RESOLUTION multiple."""
|
||
r = SNAP_RESOLUTION
|
||
return f"(round(({formula}) / {r}) * {r})"
|
||
|
||
|
||
def fetch_all_lagged_profiles(conn, contacts: list[dict]) -> dict[tuple[str, int], dict]:
|
||
"""Fetch all lagged HRRR profiles efficiently.
|
||
|
||
Strategy: Uses a temp table with all (snapped_lat, snapped_lon, lead, base_ts)
|
||
tuples, then joins against hrrr_profiles in a single query.
|
||
|
||
Returns: {(contact_id_str, lead_hours): profile_row} where profile_row is a dict
|
||
with the weather fields.
|
||
"""
|
||
if not contacts:
|
||
return {}
|
||
|
||
print(f" fetching lagged profiles for {len(contacts)} contacts at leads {LEAD_TIMES}…", file=sys.stderr)
|
||
|
||
# Build lookup rows: (contact_id, snapped_lat, snapped_lon, lead_hours, base_ts)
|
||
# We generate this in Python and pass to PostgreSQL via a VALUES clause
|
||
|
||
rows = []
|
||
for c in contacts:
|
||
slat = snap_coord(c["lat"])
|
||
slon = snap_coord(c["lon"])
|
||
base_ts = c["qso_timestamp"]
|
||
for lead in LEAD_TIMES:
|
||
rows.append((c["id"], slat, slon, lead, base_ts))
|
||
|
||
if not rows:
|
||
return {}
|
||
|
||
# Use a temp table to batch the profile lookups efficiently.
|
||
# Without ON COMMIT DROP since autocommit=True.
|
||
with conn.cursor() as cur:
|
||
cur.execute("DROP TABLE IF EXISTS _forecast_lookups")
|
||
cur.execute("""
|
||
CREATE TEMP TABLE _forecast_lookups (
|
||
contact_id text,
|
||
snapped_lat double precision,
|
||
snapped_lon double precision,
|
||
lead_hours int,
|
||
base_ts timestamptz
|
||
)
|
||
""")
|
||
# Insert in batches of 5000 rows to stay under param limits
|
||
batch_size = 5000
|
||
for i in range(0, len(rows), batch_size):
|
||
batch = rows[i:i + batch_size]
|
||
values_clauses = []
|
||
params = []
|
||
for contact_id, slat, slon, lead, base_ts in batch:
|
||
cid_str = str(contact_id)
|
||
values_clauses.append("(%s::text, %s::double precision, %s::double precision, %s::int, %s::timestamptz)")
|
||
params.extend([cid_str, slat, slon, lead, base_ts])
|
||
sql = "INSERT INTO _forecast_lookups (contact_id, snapped_lat, snapped_lon, lead_hours, base_ts) VALUES " + ", ".join(values_clauses)
|
||
cur.execute(sql, params)
|
||
|
||
# Now query hrrr_profiles joining against the lookup table
|
||
print(f" querying hrrr_profiles with {len(rows)} lookup tuples…", file=sys.stderr)
|
||
with conn.cursor(row_factory=dict_row) as cur:
|
||
cur.execute("""
|
||
SELECT
|
||
fl.contact_id,
|
||
fl.lead_hours,
|
||
h.surface_temp_c,
|
||
h.surface_dewpoint_c,
|
||
h.surface_pressure_mb,
|
||
h.pwat_mm,
|
||
h.hpbl_m,
|
||
h.min_refractivity_gradient
|
||
FROM _forecast_lookups fl
|
||
JOIN hrrr_profiles h
|
||
ON h.lat = fl.snapped_lat
|
||
AND h.lon = fl.snapped_lon
|
||
AND h.valid_time = date_trunc('hour', fl.base_ts) - make_interval(hours => fl.lead_hours)
|
||
""")
|
||
result_rows = cur.fetchall()
|
||
|
||
# Cleanup
|
||
with conn.cursor() as cur:
|
||
cur.execute("DROP TABLE IF EXISTS _forecast_lookups")
|
||
|
||
result: dict[tuple[str, int], dict] = {}
|
||
for r in result_rows:
|
||
key = (r["contact_id"], r["lead_hours"])
|
||
result[key] = {
|
||
"surface_temp_c": r["surface_temp_c"],
|
||
"surface_dewpoint_c": r["surface_dewpoint_c"],
|
||
"surface_pressure_mb": r["surface_pressure_mb"],
|
||
"pwat_mm": r["pwat_mm"],
|
||
"hpbl_m": r["hpbl_m"],
|
||
"min_refractivity_gradient": r["min_refractivity_gradient"],
|
||
}
|
||
|
||
print(f" fetched {len(result)} profile rows", file=sys.stderr)
|
||
|
||
# Report coverage per lead
|
||
for lead in LEAD_TIMES:
|
||
found = sum(1 for (cid, lh) in result if lh == lead)
|
||
print(f" lead-{lead}h: {found}/{len(contacts)} profiles found", file=sys.stderr)
|
||
|
||
return result
|
||
|
||
|
||
# ── Main ───────────────────────────────────────────────────────────────────────
|
||
|
||
def main() -> int:
|
||
dsn = os.environ.get("PROP_PROD_DB_URL")
|
||
if not dsn:
|
||
print("error: PROP_PROD_DB_URL not set", file=sys.stderr)
|
||
return 2
|
||
|
||
today = dt.date.today().isoformat()
|
||
report_dir = REPORT_DIR
|
||
report_dir.mkdir(parents=True, exist_ok=True)
|
||
|
||
# ── Load weights ────────────────────────────────────────────────────────
|
||
print("• loading band weights", file=sys.stderr)
|
||
if not WEIGHTS_PATH.exists():
|
||
print(f" warning: {WEIGHTS_PATH} not found — using global defaults", file=sys.stderr)
|
||
weights_json = {"global_weights": {
|
||
"humidity": 0.1262, "time_of_day": 0.038, "td_depression": 0.101,
|
||
"refractivity": 0.0986, "sky": 0.0841, "season": 0.1134,
|
||
"wind": 0.0841, "rain": 0.1431, "pwat": 0.1147, "pressure": 0.0967,
|
||
}, "band_overrides": {}}
|
||
else:
|
||
with open(WEIGHTS_PATH) as f:
|
||
weights_json = json.load(f)
|
||
|
||
global_weights: dict[str, float] = weights_json["global_weights"]
|
||
band_overrides: dict[str, dict] = weights_json.get("band_overrides", {})
|
||
|
||
def get_weights(band_mhz: int) -> dict[str, float]:
|
||
key = str(band_mhz)
|
||
if key in band_overrides and "weights" in band_overrides[key]:
|
||
return band_overrides[key]["weights"]
|
||
return dict(global_weights)
|
||
|
||
# ── Load data ───────────────────────────────────────────────────────────
|
||
print("• connecting to database", file=sys.stderr)
|
||
with psycopg.connect(dsn, autocommit=True) as conn:
|
||
with conn.cursor() as cur:
|
||
cur.execute("SET statement_timeout = '30min'")
|
||
contacts = load_contacts(conn)
|
||
|
||
if not contacts:
|
||
print("error: no contacts loaded — check DB connection and data", file=sys.stderr)
|
||
return 3
|
||
|
||
print(f"• loaded {len(contacts)} contacts (≥2020-01-01)", file=sys.stderr)
|
||
|
||
# ── Fetch lagged profiles ──────────────────────────────────────────
|
||
print("• fetching lagged HRRR profiles", file=sys.stderr)
|
||
all_profiles = fetch_all_lagged_profiles(conn, contacts)
|
||
|
||
# ── Score each contact at each lag ──────────────────────────────────────
|
||
print("• scoring contacts at each lead time", file=sys.stderr)
|
||
# Track which contacts have ALL required profiles (lag-0 through lag-24)
|
||
scored_by_band: dict[int, dict[int, list]] = defaultdict(lambda: defaultdict(list))
|
||
# {band_mhz: {lead_hours: [(score, distance), ...]}}
|
||
|
||
contacts_with_all = 0
|
||
for c in contacts:
|
||
band_mhz = c["band"]
|
||
band_cfg = BAND_CONFIGS.get(band_mhz)
|
||
if band_cfg is None:
|
||
continue
|
||
|
||
w = get_weights(band_mhz)
|
||
cid_str = str(c["id"])
|
||
|
||
# Check if all profiles exist
|
||
has_all = True
|
||
scores_at_lag = {}
|
||
for lead in LEAD_TIMES:
|
||
key = (cid_str, lead)
|
||
profile = all_profiles.get(key)
|
||
if profile is None:
|
||
has_all = False
|
||
break
|
||
score = compute_composite_score(profile, c, band_cfg, w)
|
||
if score is None:
|
||
has_all = False
|
||
break
|
||
scores_at_lag[lead] = (score, c["dist"])
|
||
|
||
if has_all:
|
||
contacts_with_all += 1
|
||
for lead, (score, dist) in scores_at_lag.items():
|
||
scored_by_band[band_mhz][lead].append((score, dist))
|
||
|
||
print(f"• contacts with all {len(LEAD_TIMES)} profiles: {contacts_with_all}/{len(contacts)}", file=sys.stderr)
|
||
|
||
# ── Per-band per-lag Spearman ρ ────────────────────────────────────────
|
||
print("• computing per-band per-lag Spearman ρ", file=sys.stderr)
|
||
per_band_results = {}
|
||
for band_mhz in sorted(scored_by_band.keys()):
|
||
by_lead = scored_by_band[band_mhz]
|
||
# All leads must have the same N (contacts with all profiles)
|
||
ns = [len(pairs) for pairs in by_lead.values()]
|
||
n = ns[0] if ns else 0
|
||
|
||
if n < 30:
|
||
print(f" band {band_mhz} MHz: {n} contacts — skipping (<30 threshold)", file=sys.stderr)
|
||
continue
|
||
if not all(x == n for x in ns):
|
||
print(f" band {band_mhz} MHz: inconsistent counts — skipping", file=sys.stderr)
|
||
continue
|
||
|
||
rhos = {}
|
||
for lead in LEAD_TIMES:
|
||
pairs = by_lead.get(lead, [])
|
||
scores = [s for s, d in pairs]
|
||
dists = [d for s, d in pairs]
|
||
rhos[lead] = spearman_rho(scores, dists)
|
||
|
||
rho_0 = rhos.get(0, 0.0)
|
||
per_band_results[band_mhz] = {
|
||
"n": n,
|
||
"rhos": rhos,
|
||
"delta_6h": round(rho_0 - rhos.get(6, rho_0), 4),
|
||
"delta_24h": round(rho_0 - rhos.get(24, rho_0), 4),
|
||
}
|
||
print(f" band {band_mhz} MHz: n={n}, ρ(0h)={rhos.get(0, 0):+.4f}, "
|
||
f"ρ(6h)={rhos.get(6, 0):+.4f}, ρ(24h)={rhos.get(24, 0):+.4f}", file=sys.stderr)
|
||
|
||
# ── Build report ────────────────────────────────────────────────────────
|
||
generated_at = dt.datetime.now(dt.timezone.utc).replace(microsecond=0)
|
||
|
||
json_payload = {
|
||
"generated_at": generated_at.isoformat(timespec="seconds"),
|
||
"generated_by": "scripts/validate_forecast.py",
|
||
"schema_version": 1,
|
||
"summary": {
|
||
"total_contacts": len(contacts),
|
||
"contacts_with_all_profiles": contacts_with_all,
|
||
"lead_times_hours": LEAD_TIMES,
|
||
"bands_evaluated": len(per_band_results),
|
||
},
|
||
"per_band": {
|
||
str(b): {
|
||
"n": info["n"],
|
||
"rhos": {str(k): round(v, 4) for k, v in info["rhos"].items()},
|
||
"delta_6h": info["delta_6h"],
|
||
"delta_24h": info["delta_24h"],
|
||
}
|
||
for b, info in sorted(per_band_results.items())
|
||
},
|
||
}
|
||
|
||
json_path = report_dir / f"forecast-{today}.json"
|
||
json_path.write_text(json.dumps(json_payload, indent=2, sort_keys=False) + "\n")
|
||
print(f"• wrote {json_path}", file=sys.stderr)
|
||
|
||
# ── Render Markdown table ───────────────────────────────────────────────
|
||
md = []
|
||
md.append(f"# Forecast Skill Degradation — {today}\n")
|
||
md.append(
|
||
"> Auto-generated by `scripts/validate_forecast.py`. Measures how well "
|
||
"the propagation scoring algorithm predicts contact distances as the "
|
||
"weather profile ages. The lag-0h profile uses actual conditions at "
|
||
"contact time (ground truth). Lag-Nh uses the weather profile from "
|
||
"N hours earlier — representing the forecast that would have been "
|
||
"available at that lead time.\n"
|
||
)
|
||
md.append(
|
||
f"- **Total contacts** (≥2020): {len(contacts):,}"
|
||
)
|
||
md.append(
|
||
f"- **Contacts with all {len(LEAD_TIMES)} profiles**: {contacts_with_all:,} " +
|
||
f"({contacts_with_all * 100 // max(1, len(contacts))}%)"
|
||
)
|
||
md.append(f"- **Generated**: {generated_at.isoformat(timespec='seconds')}\n")
|
||
|
||
md.append(
|
||
"Higher ρ = better distance prediction. Δρ(0h→Nh) = ρ(0h) − ρ(Nh) "
|
||
"measures skill degradation — positive values mean the forecast skill "
|
||
"is worse than using current conditions.\n"
|
||
)
|
||
|
||
md.append("## Per-Band Per-Lag Spearman ρ\n")
|
||
|
||
# Build header
|
||
col_headers = "| Band | N |"
|
||
col_sep = "|------|--:|"
|
||
for lead in LEAD_TIMES:
|
||
col_headers += f" ρ({lead}h) |"
|
||
col_sep += "-------:|" if lead == 0 else "------:|"
|
||
col_headers += " Δρ(0→6h) | Δρ(0→24h) |"
|
||
col_sep += "----------|------------|"
|
||
|
||
md.append(col_headers)
|
||
md.append(col_sep)
|
||
|
||
# Sort bands by contact count descending for display
|
||
sorted_bands = sorted(per_band_results.keys(), key=lambda b: per_band_results[b]["n"], reverse=True)
|
||
for band_mhz in sorted_bands:
|
||
info = per_band_results[band_mhz]
|
||
row = f"| {band_mhz} MHz | {info['n']} |"
|
||
for lead in LEAD_TIMES:
|
||
rho = info["rhos"].get(lead, 0.0)
|
||
row += f" {rho:+.4f} |"
|
||
row += f" {info['delta_6h']:+.4f} | {info['delta_24h']:+.4f} |"
|
||
md.append(row)
|
||
|
||
# ── Key findings ────────────────────────────────────────────────────────
|
||
md.append("")
|
||
md.append("## Key Findings\n")
|
||
|
||
# Check monotonic degradation
|
||
monotonic_bands = 0
|
||
non_monotonic_bands = []
|
||
for band_mhz in sorted_bands:
|
||
info = per_band_results[band_mhz]
|
||
rhos = info["rhos"]
|
||
prev = None
|
||
monotonic = True
|
||
for lead in LEAD_TIMES:
|
||
rho = rhos.get(lead, 0.0)
|
||
if prev is not None and rho > prev + 0.001: # slight tolerance
|
||
monotonic = False
|
||
prev = rho
|
||
if monotonic:
|
||
monotonic_bands += 1
|
||
else:
|
||
non_monotonic_bands.append(band_mhz)
|
||
|
||
md.append(f"- **Monotonic degradation** (ρ decreases with lead time): {monotonic_bands}/{len(sorted_bands)} bands")
|
||
|
||
if non_monotonic_bands:
|
||
md.append(f"- **Non-monotonic bands** (ρ increases at some lead): {', '.join(f'{b} MHz' for b in non_monotonic_bands)}")
|
||
md.append(
|
||
" - This may indicate the algorithm puts disproportionate weight "
|
||
"on noisy short-term features that vary within hours. If ρ(6h) > ρ(0h), "
|
||
"the current-conditions score is noisier than the 6h-lagged score — "
|
||
"weather features at lag-0 may introduce variance that degrades the "
|
||
"correlation."
|
||
)
|
||
|
||
# @6h vs @0h comparison
|
||
rho_6h_better = 0
|
||
for band_mhz in sorted_bands:
|
||
info = per_band_results[band_mhz]
|
||
r0 = info["rhos"].get(0, 0.0)
|
||
r6 = info["rhos"].get(6, 0.0)
|
||
if r6 > r0:
|
||
rho_6h_better += 1
|
||
if rho_6h_better > 0:
|
||
md.append(f"- **ρ(6h) > ρ(0h)**: {rho_6h_better} bands — the 6-hour lagged score outperforms current-conditions")
|
||
|
||
# Average degradation
|
||
if per_band_results:
|
||
avg_delta_6 = sum(info["delta_6h"] for info in per_band_results.values()) / len(per_band_results)
|
||
avg_delta_24 = sum(info["delta_24h"] for info in per_band_results.values()) / len(per_band_results)
|
||
md.append(f"- **Mean Δρ(0→6h)**: {avg_delta_6:+.4f}")
|
||
md.append(f"- **Mean Δρ(0→24h)**: {avg_delta_24:+.4f}")
|
||
else:
|
||
md.append("- **No bands passed the ≥30 contact threshold.**")
|
||
|
||
# Caveats
|
||
md.append("")
|
||
md.append("## Caveats\n")
|
||
md.append(
|
||
"- **Lagged profiles use analysis data**: The HRRR profiles at T−N "
|
||
"represent actual weather conditions at that earlier time, not a true "
|
||
"forecast run initialized at T−N. Real HRRR forecasts would contain "
|
||
"model error growth; these measurements capture the degradation from "
|
||
"weather *evolution* alone (the lower bound of forecast skill loss)."
|
||
)
|
||
md.append(
|
||
"- **Snap to nearest grid point**: Profiles are looked up at the HRRR "
|
||
"grid point nearest to the contact location (rounded to 0.125°). "
|
||
"A single point approximates the full path's conditions."
|
||
)
|
||
md.append(
|
||
"- **Filter threshold**: Only bands with ≥30 contacts having ALL "
|
||
"6 profiles (lag-0 through lag-24) are included in the table."
|
||
)
|
||
md.append(
|
||
"- **Censored data**: Contacts are confirmed success events only. "
|
||
"Correlation ρ measures discrimination among successful contacts."
|
||
)
|
||
md.append("")
|
||
|
||
md_path = report_dir / f"forecast-{today}.md"
|
||
md_text = "\n".join(md)
|
||
md_path.write_text(md_text + "\n")
|
||
print(f"• wrote {md_path}", file=sys.stderr)
|
||
|
||
# ── Print to stdout ─────────────────────────────────────────────────────
|
||
print(md_text)
|
||
|
||
return 0
|
||
|
||
|
||
if __name__ == "__main__":
|
||
raise SystemExit(main())
|