#!/usr/bin/env python3
"""Tier-1 probe: is the 2026-06-30 Gaetano fire's thermal signal present-but-
sub-threshold in our SEVIRI/FCI data, or truly absent?

Gold truth (forest-ranger/Gaetano confirmed), Agrigento W Sicily:
  p1 = 37.4602054, 13.2264283   p2 = 37.456807, 13.224735
  fire observed ~19:20 UTC 30 Jun 2026 (21:20 CEST). ALL feeds missed it;
  sky was clear (ECMWF tcc 0.06-0.20 @18:00Z) so this is not occlusion.

What this does (single-pixel, reuses the LIVE repo detectors — no reimpl):
  For each cached SEVIRI MSG15-RSS .nat frame in the burn window, load
  bt_mir(3.9um)/bt_tir(10.8um), find the grid pixel nearest each truth point,
  and record:
    - bt_mir, bt_tir, delta at the truth pixel + 3x3 max
    - baseline z-score  (bt_mir - mean_mir)/std_mir  vs the per-pixel-hour-of-
      day diurnal baseline (<per-pixel diurnal baseline .npz>)
    - whether detect_subpixel_v1 fires at that pixel at DEFAULT thresholds
    - ... and at CITIZEN-PRIMED RELAXED thresholds (mir 290K, z 2.0) — i.e. the
      relaxation the ground-report prior now applies at this exact spot.

Verdict:
  * a warm MIR anomaly (z rising >~2, or bt_mir climbing above the diurnal
    baseline) at the truth pixel during the burn window  => SIGNAL PRESENT:
    tune the floor / the citizen prior would catch it.
  * flat-at-baseline the whole window  => ABSENT below SEVIRI's noise floor:
    thermal geostationary cannot get this fire; DNB / next-day S2 SWIR are
    mandatory. (Run the DNB half separately on the DOY182 granule.)

Venue-agnostic: runs in the repo env (imports src.*). Footprint is a few
minutes over ~40 frames, so it runs in a few minutes in any batch environment.

Usage:
  python -m probe_geostationary_floor_seviri \
      --seviri-cache $EUMETSAT_CACHE/seviri_rss_cache \
      --start 2026-06-30T19:00 --end 2026-07-01T06:00 \
      --out $EUMETSAT_CACHE/probe_gaetano_20260630.json
"""
from __future__ import annotations

import argparse
import glob
import json
import logging
import os
import re
from datetime import datetime, timezone
from pathlib import Path

import numpy as np

logging.basicConfig(level=logging.INFO,
                    format="%(asctime)s [%(levelname)s] %(name)s: %(message)s")
log = logging.getLogger("probe_gaetano")

TRUTH = [("p1", 37.4602054, 13.2264283), ("p2", 37.456807, 13.224735)]
# SEVIRI RSS .nat filename carries the sensing timestamp: ...-<YYYYMMDDHHMMSS>...Z
_TS_RE = re.compile(r"-(\d{14})\.\d+Z-")


def _fname_ts(path: str):
    m = _TS_RE.search(os.path.basename(path))
    if not m:
        return None
    return datetime.strptime(m.group(1), "%Y%m%d%H%M%S").replace(tzinfo=timezone.utc)


def _nearest_idx(lat_grid, lon_grid, lat0, lon0):
    d2 = (lat_grid - lat0) ** 2 + ((lon_grid - lon0) * np.cos(np.radians(lat0))) ** 2
    return np.unravel_index(int(np.nanargmin(d2)), d2.shape)


def _win3x3(arr, r, c):
    r0, r1 = max(0, r - 1), min(arr.shape[0], r + 2)
    c0, c1 = max(0, c - 1), min(arr.shape[1], c + 2)
    return arr[r0:r1, c0:c1]


def main(argv=None) -> int:
    ap = argparse.ArgumentParser(prog="probe_gaetano_fire")
    ap.add_argument("--seviri-cache", required=True)
    ap.add_argument("--start", default="2026-06-30T19:00")
    ap.add_argument("--end", default="2026-07-01T06:00")
    ap.add_argument("--out", required=True)
    args = ap.parse_args(argv)

    start = datetime.fromisoformat(args.start).replace(tzinfo=timezone.utc)
    end = datetime.fromisoformat(args.end).replace(tzinfo=timezone.utc)

    # Reuse the LIVE repo code: reader, detector, per-pixel diurnal baseline.
    from src.data_sources.seviri_rss import SeviriRssClient
    from src.detectors import subpixel_v1 as sp

    client = SeviriRssClient()
    baseline = None  # lazy-loaded once we know the grid shape

    nats = sorted(glob.glob(os.path.join(args.seviri_cache, "MSG4-SEVI-MSG15-*.nat")))
    frames = [(p, _fname_ts(p)) for p in nats]
    frames = [(p, t) for (p, t) in frames if t is not None and start <= t <= end]
    log.info("probe: %d SEVIRI frames in window %s .. %s", len(frames), start, end)

    series = []
    for path, ts in frames:
        try:
            fr = client._read_and_crop(Path(path))
        except Exception as exc:  # noqa: BLE001
            log.warning("read failed %s: %s", os.path.basename(path), exc)
            continue
        if baseline is None:
            baseline = sp._load_baseline(shape=fr.bt_mir.shape)

        rec = {"ts": (fr.scan_start or ts).isoformat(),
               "product_id": fr.product_id, "points": {}}

        # Default-threshold and relaxed (citizen-primed) detections for this frame.
        det_default = sp.detect_subpixel_v1(
            fr.bt_mir, fr.bt_tir, fr.lat, fr.lon,
            source_stream="seviri_rss", timestamp=fr.scan_start or ts,
            baseline=baseline)
        det_relaxed = sp.detect_subpixel_v1(
            fr.bt_mir, fr.bt_tir, fr.lat, fr.lon,
            source_stream="seviri_rss", timestamp=fr.scan_start or ts,
            mir_threshold_k=290.0, lightning_primed_threshold_k=288.0,
            baseline_z_threshold=2.0, baseline=baseline)

        for name, la, lo in TRUTH:
            r, c = _nearest_idx(fr.lat, fr.lon, la, lo)
            bt_mir = float(fr.bt_mir[r, c])
            bt_tir = float(fr.bt_tir[r, c])
            z = None
            if baseline is not None and "mean_mir" in baseline and "std_mir" in baseline:
                try:
                    hod = (fr.scan_start or ts).hour
                    mean = baseline["mean_mir"]
                    std = baseline["std_mir"]
                    m = mean[hod, r, c] if mean.ndim == 3 else mean[r, c]
                    s = std[hod, r, c] if std.ndim == 3 else std[r, c]
                    if s and np.isfinite(s) and s > 0:
                        z = (bt_mir - float(m)) / float(s)
                except Exception:
                    z = None

            def _hit(dets):
                for d in dets:
                    if abs(getattr(d, "lat", 1e9) - la) < 0.03 and \
                       abs(getattr(d, "lng", getattr(d, "lon", 1e9)) - lo) < 0.03:
                        return True
                return False

            rec["points"][name] = {
                "bt_mir_k": round(bt_mir, 3),
                "bt_tir_k": round(bt_tir, 3),
                "delta_k": round(bt_mir - bt_tir, 3),
                "bt_mir_3x3_max_k": round(float(np.nanmax(_win3x3(fr.bt_mir, r, c))), 3),
                "baseline_z": round(z, 3) if z is not None else None,
                "subpixel_default_fired": _hit(det_default),
                "subpixel_relaxed_fired": _hit(det_relaxed),
            }
        series.append(rec)
        log.info("%s  p1 bt_mir=%.2f z=%s", rec["ts"],
                 rec["points"]["p1"]["bt_mir_k"], rec["points"]["p1"]["baseline_z"])

    # Verdict heuristic over the burn window.
    def _mx(name, key):
        vals = [s["points"][name][key] for s in series
                if s["points"][name].get(key) is not None]
        return max(vals) if vals else None

    verdict = {}
    for name, _, _ in TRUTH:
        zmax = _mx(name, "baseline_z")
        dmax = _mx(name, "delta_k")
        any_relaxed = any(s["points"][name]["subpixel_relaxed_fired"] for s in series)
        any_default = any(s["points"][name]["subpixel_default_fired"] for s in series)
        present = (zmax is not None and zmax >= 2.0) or any_relaxed or any_default
        verdict[name] = {
            "max_baseline_z": zmax, "max_delta_k": dmax,
            "relaxed_fired_any_frame": any_relaxed,
            "default_fired_any_frame": any_default,
            "signal_present": bool(present),
            "call": ("PRESENT-SUBTHRESHOLD (tune floor / citizen prior catches it)"
                     if present else
                     "ABSENT below SEVIRI floor (DNB / next-day S2 SWIR mandatory)"),
        }

    out = {"truth_points": TRUTH, "window": [args.start, args.end],
           "n_frames": len(series), "verdict": verdict, "series": series}
    Path(args.out).write_text(json.dumps(out, indent=2, default=str))
    log.info("VERDICT: %s", json.dumps(verdict, indent=2, default=str))
    log.info("wrote %s", args.out)
    return 0


if __name__ == "__main__":
    raise SystemExit(main())