#!/usr/bin/env python3 """ Train LightGBM models for HK weather prediction targets. The training data simulates the relationship between NWP model forecasts and actual observations. NWP models have systematic errors: - Temperature: RMSE ~1.5°C at 24h lead - Precipitation probability: poor calibration, often overconfident - Wind: RMSE 3-5 km/h at 24h lead The model learns to MAP noisy forecast features → binary outcome truth. """ import argparse import sys import json from datetime import datetime, timedelta from pathlib import Path from typing import Optional, Dict, Tuple import numpy as np import pandas as pd sys.path.insert(0, str(Path(__file__).parent.parent)) from ml.features import FeatureEngine from ml.model import WeatherModel, TARGET_DEFINITIONS, MODEL_DIR from config import HK_COORDS def _add_nwp_forecast_error( daily_truth: pd.DataFrame, hourly_truth: pd.DataFrame, rng: np.random.RandomState, ) -> Tuple[pd.DataFrame, pd.DataFrame]: """ Add realistic NWP forecast errors to truth data. Returns (daily_forecast, hourly_forecast) simulating: - Temperature: RMSE 1.5-2.5°C, warm bias in summer anticyclones - Precipitation: continuous calibrated probabilities (not just 0/100), systematic overforecasting of light rain, underforecasting of heavy - Wind: multiplicative errors 0.7-1.5x - Cloud cover: RMSE 15-20% """ n_days = len(daily_truth) # === TEMPERATURE: larger noise for wider training distribution === temp_error_max = rng.normal(0.3, 2.0, n_days) # μ=0.3 bias, σ=2.0 RMSE temp_error_min = rng.normal(0.2, 1.8, n_days) # Random injection of larger errors (10% of days have outlier errors) outlier_mask = rng.random(n_days) < 0.10 temp_error_max[outlier_mask] += rng.normal(0, 3.0, outlier_mask.sum()) temp_error_min[outlier_mask] += rng.normal(0, 2.5, outlier_mask.sum()) daily_fc = daily_truth.copy() if "temperature_2m_max" in daily_fc.columns: daily_fc["temperature_2m_max"] = np.clip(daily_truth["temperature_2m_max"] + temp_error_max, 5, 42) if "temperature_2m_min" in daily_fc.columns: daily_fc["temperature_2m_min"] = np.clip(daily_truth["temperature_2m_min"] + temp_error_min, 0, 33) # === PRECIPITATION: continuous calibrated probabilities + multiplicative rain error === # NWP models output continuous probabilities, not 0/100 # Base prob from truth, then add calibration noise true_prob = daily_truth["precipitation_probability_max"].values / 100.0 # Systematic miscalibration: NWP overestimates low prob, underestimates high prob calibration_bias = 0.15 * (0.5 - true_prob) # +7.5% at prob=0, -7.5% at prob=1 calibration_noise = rng.normal(0, 0.15, n_days) fc_prob = np.clip(true_prob + calibration_bias + calibration_noise, 0.01, 0.99) if "precipitation_probability_max" in daily_fc.columns: daily_fc["precipitation_probability_max"] = fc_prob * 100.0 # Rain amount: multiplicative error, more noise on heavy rain rain_mult_error = np.where( daily_truth["precipitation_sum"] > 5, rng.lognormal(0, 0.4, n_days), # High variance for heavy rain rng.lognormal(0, 0.25, n_days), # Lower variance for light rain ) if "precipitation_sum" in daily_fc.columns: daily_fc["precipitation_sum"] = daily_truth["precipitation_sum"] * rain_mult_error # === WIND: multiplicative with 10% outlier days === wind_mult = rng.lognormal(0, 0.20, n_days) gust_mult = rng.lognormal(0, 0.30, n_days) outlier_wind = rng.random(n_days) < 0.10 wind_mult[outlier_wind] *= rng.uniform(1.3, 2.0, outlier_wind.sum()) gust_mult[outlier_wind] *= rng.uniform(1.3, 2.5, outlier_wind.sum()) for col, mult in [("wind_speed_10m_max", wind_mult), ("wind_gusts_10m_max", gust_mult)]: if col in daily_fc.columns: daily_fc[col] = np.clip(daily_truth[col] * mult, 0, 200) # === CLOUD COVER: systematic bias (underestimate in convective conditions) === if "weather_code" in daily_fc.columns: daily_fc["weather_code"] = daily_truth["weather_code"] # === HOURLY: larger noise ranges === hourly_fc = hourly_truth.copy() n_hours = len(hourly_fc) # Temperature: diurnal-cycle-aware errors (larger at night) hour_of_day = np.array([i % 24 for i in range(n_hours)]) t_noise_scale = 1.2 + 0.8 * np.sin(2 * np.pi * (hour_of_day - 14) / 24) # Peak error at night if "temperature_2m" in hourly_fc.columns: hourly_fc["temperature_2m"] = np.clip( hourly_truth["temperature_2m"] + rng.normal(0, 2.0, n_hours) * t_noise_scale, -5, 45, ) # Humidity: large errors (NWP struggles with boundary layer moisture) if "relative_humidity_2m" in hourly_fc.columns: rh_err = rng.normal(-3, 12, n_hours) # More error during convective hours convective_mask = (hour_of_day > 11) & (hour_of_day < 19) rh_err[convective_mask] *= 1.5 hourly_fc["relative_humidity_2m"] = np.clip(hourly_truth["relative_humidity_2m"] + rh_err, 15, 100) # Cloud cover: large RMSE if "cloud_cover" in hourly_fc.columns: hourly_fc["cloud_cover"] = np.clip(hourly_truth["cloud_cover"] + rng.normal(0, 20, n_hours), 0, 100) for level in ["cloud_cover_low", "cloud_cover_mid", "cloud_cover_high"]: if level in hourly_fc.columns: hourly_fc[level] = np.clip(hourly_truth[level] + rng.normal(0, 15, n_hours), 0, 100) # Pressure: typical errors if "surface_pressure" in hourly_fc.columns: hourly_fc["surface_pressure"] = hourly_truth["surface_pressure"] + rng.normal(0, 3.0, n_hours) # Wind: multiplicative for col in ["wind_speed_10m", "wind_speed_100m", "wind_gusts_10m"]: if col in hourly_fc.columns: mult = rng.lognormal(0, 0.25, n_hours) hourly_fc[col] = hourly_truth[col] * mult return daily_fc, hourly_fc def bootstrap_training_data(n_samples: int = 8000) -> Tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame, pd.DataFrame]: """ Generate NWP forecast + observation training pairs. Returns (daily_forecast, hourly_forecast, daily_truth, hourly_truth) where forecast has realistic NWP errors and truth is the actual observation. """ rng = np.random.RandomState(42) rng_noise = np.random.RandomState(99) start_date = datetime(2015, 1, 1) dates = [start_date + timedelta(days=i) for i in range(n_samples)] doy = np.array([d.timetuple().tm_yday for d in dates]) # === Generate OBSERVATION TRUTH (clean, no NWP error) === tmax_base = 26.0 + 7.0 * np.sin(2 * np.pi * (doy - 200) / 365.25) tmax = np.clip(tmax_base + rng_noise.normal(0, 2.0, n_samples), 8, 38) tmin = np.clip(tmax - (7.0 + rng_noise.exponential(2.0, n_samples)), 4, 30) tmean = (tmax + tmin) / 2 rain_seasonal = 4.0 + 10.0 * np.maximum(0, np.sin(2 * np.pi * (doy - 172) / 365.25)) rain_day_mask_prob = 0.3 + 0.4 * np.maximum(0, np.sin(2 * np.pi * (doy - 172) / 365.25)) rain_day = rng_noise.random(n_samples) < rain_day_mask_prob rain_sum = np.where(rain_day, rng_noise.exponential(rain_seasonal, n_samples), 0) rain_sum = np.where(rain_sum < 0.1, 0, rain_sum) # Continuous precipitation probability: beta distribution centered on actual prob from scipy.stats import beta as beta_dist precip_prob = np.zeros(n_samples) for i in range(n_samples): # Center the beta around the climatological rain prob p = rain_day_mask_prob[i] a = max(0.5, p * 8) b = max(0.5, (1 - p) * 8) precip_prob[i] = rng_noise.beta(a, b) * 100.0 precip_prob = np.clip(precip_prob, 0.5, 99.5) wind_base = 15 + 8 * np.maximum(0, np.sin(2 * np.pi * (doy - 180) / 365.25)) wind_max = np.clip(wind_base + rng_noise.exponential(5, n_samples), 3, 120) gusts_max = np.clip(wind_max * (1.0 + rng_noise.exponential(0.5, n_samples)), wind_max, 200) wind_dir = rng_noise.uniform(0, 360, n_samples) cloud = np.clip(30 + rng_noise.beta(2, 3, n_samples) * 70 * (0.5 + 0.5 * (rain_sum > 0)), 0, 100) pressure = np.clip(1013 - 5 * np.sin(2 * np.pi * (doy - 172) / 365.25) + rng_noise.normal(0, 3, n_samples), 980, 1035) sw_rad = np.clip(5.0 + 10.0 * np.sin(2 * np.pi * (doy - 172) / 365.25) * (1 - cloud / 100) + rng_noise.normal(0, 2, n_samples), 0, 30) daily_truth = pd.DataFrame({ "temperature_2m_max": tmax, "temperature_2m_min": tmin, "temperature_2m_mean": tmean, "precipitation_sum": rain_sum, "precipitation_probability_max": precip_prob, "rain_sum": rain_sum, "wind_speed_10m_max": wind_max, "wind_gusts_10m_max": gusts_max, "wind_direction_10m_dominant": wind_dir, "shortwave_radiation_sum": sw_rad, "et0_fao_evapotranspiration": sw_rad * 0.4, "weather_code": np.where(rain_sum > 0, np.where(rain_sum > 10, 63, 61), 0), }, index=pd.to_datetime(dates)) # Hourly truth hours_per_day = 24 total_hours = n_samples * hours_per_day hour_of_day = np.tile(np.arange(24), n_samples) t_range = np.repeat(tmax - tmin, hours_per_day) t_phase = 2 * np.pi * (hour_of_day - 14) / 24 daily_rh = np.clip(78 + 12 * np.sin(2 * np.pi * (doy - 172) / 365.25) + rng_noise.normal(0, 6, n_samples), 45, 98) dewpoint = tmean - ((100 - daily_rh) / 5.0) + rng_noise.normal(0, 0.5, n_samples) t_hourly = np.repeat(tmin, hours_per_day) + t_range * (0.5 + 0.5 * np.cos(t_phase)) + rng_noise.normal(0, 0.5, total_hours) rh_hourly = np.clip(np.repeat(daily_rh, hours_per_day) - 5 * np.cos(t_phase) + rng_noise.normal(0, 3, total_hours), 20, 100) hour_timestamps = [start_date + timedelta(hours=i) for i in range(total_hours)] hourly_truth = pd.DataFrame({ "temperature_2m": t_hourly, "relative_humidity_2m": rh_hourly, "dew_point_2m": np.repeat(dewpoint, hours_per_day) + rng_noise.normal(0, 0.5, total_hours), "apparent_temperature": t_hourly + rng_noise.exponential(2.0, total_hours), "precipitation_probability": np.clip(np.repeat(precip_prob, hours_per_day) / 24 + rng_noise.normal(0, 1, total_hours), 0, 100), "precipitation": np.repeat(rain_sum, hours_per_day) / 24 * rng_noise.uniform(0.5, 1.5, total_hours), "rain": np.repeat(rain_sum, hours_per_day) / 24, "cloud_cover": np.clip(np.repeat(cloud, hours_per_day) + rng_noise.normal(0, 5, total_hours), 0, 100), "cloud_cover_low": np.clip(np.repeat(cloud * 0.6, hours_per_day), 0, 100), "cloud_cover_mid": np.clip(np.repeat(cloud * 0.3, hours_per_day), 0, 100), "cloud_cover_high": np.clip(np.repeat(cloud * 0.2, hours_per_day), 0, 100), "wind_speed_10m": np.repeat(wind_max, hours_per_day) * 0.5 * (0.5 + 0.5 * np.cos(t_phase)), "wind_speed_100m": np.repeat(wind_max, hours_per_day) * 1.3 * 0.6, "wind_gusts_10m": np.repeat(gusts_max, hours_per_day) * (0.3 + 0.7 * rng_noise.beta(2, 5, total_hours)), "wind_direction_10m": np.repeat(wind_dir, hours_per_day) + rng_noise.normal(0, 10, total_hours), "surface_pressure": np.repeat(pressure, hours_per_day) + rng_noise.normal(0, 0.5, total_hours), "visibility": np.clip(15000 - np.repeat(rain_sum, hours_per_day) * 500 + rng_noise.normal(0, 2000, total_hours), 500, 25000), }, index=pd.to_datetime(hour_timestamps)) # === Add NWP forecast errors === daily_fc, hourly_fc = _add_nwp_forecast_error(daily_truth.copy(), hourly_truth.copy(), rng) return daily_fc, hourly_fc, daily_truth, hourly_truth def train_targets(target_names=None, n_bootstrap=8000, test_split=0.2): """Train models with proper NWP forecast → observation mapping.""" if target_names is None: target_names = list(TARGET_DEFINITIONS.keys()) print(f"Training {len(target_names)} models with realistic NWP errors") print(f" Samples: {n_bootstrap} (test: {test_split:.0%})") print() daily_fc, hourly_fc, daily_truth, hourly_truth = bootstrap_training_data(n_bootstrap) engine = FeatureEngine() # Features from FORECAST (noisy NWP output) X = engine.transform(daily_fc, hourly_fc) print(f"Features: {X.shape[1]} from NWP forecast output") split_idx = int(len(daily_fc) * (1 - test_split)) X_train, X_test = X[:split_idx], X[split_idx:] daily_truth_train, daily_truth_test = daily_truth.iloc[:split_idx], daily_truth.iloc[split_idx:] results = {} for target_name in target_names: tdef = TARGET_DEFINITIONS[target_name] print(f"\n{'='*60}") print(f" {target_name} — {tdef['description']}") print(f"{'='*60}") # Targets from TRUTH (actual observation) y_train = WeatherModel.build_target(daily_truth_train, tdef["variable"], tdef["threshold"], tdef["op"]) y_test = WeatherModel.build_target(daily_truth_test, tdef["variable"], tdef["threshold"], tdef["op"]) p_yes = y_train.mean() * 100 print(f" Class balance: {p_yes:.1f}% YES / {100-p_yes:.1f}% NO") model = WeatherModel(target_name, mode="lr") model.train(X_train, y_train, X_test, y_test) metrics = model.evaluate(X_test, y_test) model.save() results[target_name] = metrics print(f" Pre-calibration Brier: — ") print(f" Post-calibration:") print(f" Brier: {metrics['brier_score']:.4f} AUC: {metrics['roc_auc']:.3f}") print(f" Predicted mean: {metrics['p_yes_predicted']:.1f}% Actual: {metrics['p_yes_actual']:.1f}%") print(f" Calibration: {metrics['calibration_method']}") print(f" Top 8 features:") for feat, imp in list(model.top_features(8).items()): print(f" {feat:30s} {imp:>10.1f}") print(f"\n{'='*60}") print("TRAINING SUMMARY") print(f"{'='*60}") print(f"{'Target':<25s} {'Brier':>8s} {'AUC':>8s} {'Cal Err%':>9s} {'Cal':>10s}") print("-" * 65) for name, m in results.items(): cal_err = abs(m["p_yes_predicted"] - m["p_yes_actual"]) print(f"{name:<25s} {m['brier_score']:>8.4f} {m['roc_auc']:>8.3f} {cal_err:>9.1f} {m['calibration_method']:>10s}") print(f"\nModels saved to: {MODEL_DIR}") return results def main(): parser = argparse.ArgumentParser(description="Train HK weather prediction models") parser.add_argument("--target", type=str, default=None) parser.add_argument("--samples", type=int, default=8000) parser.add_argument("--all", action="store_true", default=False) args = parser.parse_args() targets = [args.target] if args.target else list(TARGET_DEFINITIONS.keys()) results = train_targets(targets, n_bootstrap=args.samples) MODEL_DIR.mkdir(parents=True, exist_ok=True) with open(MODEL_DIR / "training_summary.json", "w") as f: json.dump({ "training_date": datetime.now().isoformat(), "n_samples": args.samples, "results": results, }, f, indent=2, default=str) if __name__ == "__main__": main()