""" PHOENIX research — does the FULL raw SEVIRI spectral vector separate REAL from FALSE fire candidates, where DERIVED features (entry 0037) could not? Decisive extension of entry 0037 (derived features = temperature/FRP/persistence/confidence did NOT separate real vs false candidates). Here we test the FULL raw 11-band SEVIRI spectral vector (all narrow channels + pairwise differences) at each candidate pixel, on the large 2,591-frame corpus (entry 0058's data), with honest spatial cross-validation and a placebo. Pipeline: 1. Load all SEVIRI frames (11 narrow channels per pixel). 2. Front-end: flag candidate hot pixels by contextual MIR-TIR spatial anomaly (a Giglio/ Roberts-style contextual detector) — yields real fires + spectrally-fire-like false alarms. 3. Label each candidate by NASA FIRMS ground truth (<=0.05 deg, <=+/-20 min). 4. Three feature sets per candidate: DERIVED (0037 baseline): MIR, TIR, MIR-TIR, spatial-z, temporal-diff, local bg, etc. SPECTRAL (the new test): 11 raw band BTs + all 55 pairwise band differences. COMBINED. 5. HistGB classifier, GroupKFold spatial CV (0.5 deg blocks). AUC + precision@recall. 6. Placebo (shuffled labels) -> must be ~0.5 (no leakage). Numbers only. RunPod (EU-RO-1), phoenix-training volume. No DGX. """ import numpy as np, glob, os, re, sqlite3, time, json from datetime import datetime, timezone from sklearn.ensemble import HistGradientBoostingClassifier from sklearn.model_selection import GroupKFold from sklearn.metrics import roc_auc_score FRAMEDIR="/workspace/seviri_frames"; DB="/workspace/data/ground_truth.sqlite" OUT="/workspace/spec12_candidate_precision_results.json" FEAT_OUT="/workspace/spec12_candidate_features.npz" def log(m): print("[spec12 %s] %s"%(time.strftime("%H:%M:%S"),m),flush=True) CH=['bt_ir039','bt_wv062','bt_wv073','bt_ir087','bt_ir097','bt_ir108','bt_ir120','bt_ir134', 'bt_vis006','bt_vis008','bt_nir016'] # 11 SEVIRI narrow channels (HRV not stored) MIR='bt_ir039'; TIR='bt_ir108' FIRMS={'firms_viirs_noaa20','firms_viirs_noaa21','firms_viirs_snpp','firms_modis_nrt', 'slstr_frp_s3a','slstr_frp_s3b','mtg_af_l2','viirs_dnb_fire', 'firms_viirs_noaa20_nrt_archive','firms_viirs_noaa21_nrt_archive','firms_viirs_snpp_nrt_archive', 'firms_modis_nrt_archive','firms_viirs_noaa20_global'} FN_TS=re.compile(r'-NA-(\d{14})') def fname_ts(n): m=FN_TS.search(n); return datetime.strptime(m.group(1),'%Y%m%d%H%M%S').replace(tzinfo=timezone.utc) if m else None def parse_ts(ts): if not ts: return None s=ts.strip() for f in ('%Y-%m-%dT%H:%M:%S','%a, %d %b %Y %H:%M:%S'): try: base=s.replace('Z','').split('.')[0].split('+')[0].strip() return datetime.strptime(base,f).replace(tzinfo=timezone.utc) except: pass try: return datetime.fromisoformat(s.replace('Z','+00:00').split('.')[0]).replace(tzinfo=timezone.utc) except: return None # ---- load frames ---- files=sorted(glob.glob(os.path.join(FRAMEDIR,'*.npz'))) frames=[];ftime=[];flat=[];flon=[] for p in files: t=fname_ts(os.path.basename(p)) if t is None: continue z=np.load(p,allow_pickle=True) if MIR not in z or TIR not in z: continue try: cube=np.stack([np.nan_to_num(z[c].astype(np.float32),nan=0.0) for c in CH],0) except: continue if cube.shape[1:]!=(57,100): continue frames.append(cube);ftime.append(t);flat.append(z['lat'].astype(np.float32));flon.append(z['lon'].astype(np.float32)) F=np.stack(frames,0);LAT=np.stack(flat,0);LON=np.stack(flon,0);T,C,H,W=F.shape order=np.argsort([t.timestamp() for t in ftime]);F=F[order];LAT=LAT[order];LON=LON[order] ftime=[ftime[i] for i in order] log("frames=%d cube=%s %s..%s"%(T,(C,H,W),ftime[0].date(),ftime[-1].date())) mir=CH.index(MIR); tir=CH.index(TIR) # ---- FIRMS truth -> per (frame,pixel) ---- con=sqlite3.connect(DB);cur=con.cursor() fires=[] for (s,la,lo,ts,fr) in cur.execute("SELECT source,lat,lng,ts,frp_w FROM external_fires"): if s not in FIRMS: continue t=parse_ts(ts) if t is None or la is None or lo is None: continue fires.append((la,lo,t.timestamp(),fr or 0.0)) log("FIRMS truth rows=%d"%len(fires)) ftime_arr=np.array([t.timestamp() for t in ftime]) Y=np.zeros((T,H,W),np.uint8); FRP=np.zeros((T,H,W),np.float32) for (la,lo,tt,fr) in fires: dt=np.abs(ftime_arr-tt); fi=int(dt.argmin()) if dt[fi]>20*60: continue d2=(LAT[fi]-la)**2+(LON[fi]-lo)**2; idx=int(d2.argmin()); hi,wi=np.unravel_index(idx,(H,W)) if d2[hi,wi]>(0.05**2): continue Y[fi,hi,wi]=1; FRP[fi,hi,wi]=max(FRP[fi,hi,wi],fr) log("FIRMS-matched positive pixels=%d"%int(Y.sum())) # ---- front-end: contextual MIR-TIR spatial anomaly ---- def spatial_z(x): # per-frame z over valid (>0) pixels out=np.zeros_like(x) for i in range(x.shape[0]): v=x[i]; m=v[v>0] mu=m.mean() if m.size else 0.0; sd=m.std()+1e-6 out[i]=(v-mu)/sd return out MIRf=F[:,mir]; TIRf=F[:,tir] doz=MIRf-TIRf doz_z=spatial_z(doz); mir_z=spatial_z(MIRf) # temporal diff (frame - prev), MIR dprev=np.zeros_like(MIRf); dprev[1:]=MIRf[1:]-MIRf[:-1]; tdiff_z=spatial_z(dprev) # local background (3x3 exclude center) mean/std of MIR via shifts def local_bg(x): s=np.zeros_like(x); c=np.zeros_like(x); sq=np.zeros_like(x) for dh in (-1,0,1): for dw in (-1,0,1): if dh==0 and dw==0: continue sh=np.roll(np.roll(x,dh,axis=1),dw,axis=2) s+=sh; sq+=sh*sh; c+=1.0 mu=s/c; var=np.clip(sq/c-mu*mu,0,None); return mu,np.sqrt(var)+1e-6 bg_mu,bg_sd=local_bg(MIRf) THRESH=2.5 cand=(doz_z>=THRESH)&(MIRf>280.0) # hot contextual anomaly, physical MIR floor ci=np.where(cand) nC=len(ci[0]); log("candidates=%d (thresh doz_z>=%.1f & MIR>280K)"%(nC,THRESH)) labels=Y[ci]; log("candidate real=%d false=%d real_frac=%.4f"%(int(labels.sum()),int((1-labels).sum()),labels.mean())) # ---- features per candidate ---- ti,hi,wi=ci raw=np.stack([F[ti,k,hi,wi] for k in range(C)],1) # (nC,11) raw band BTs # pairwise differences pairs=[(a,b) for a in range(C) for b in range(a+1,C)] diffs=np.stack([raw[:,a]-raw[:,b] for (a,b) in pairs],1) # (nC,55) SPECTRAL=np.concatenate([raw,diffs],1) # 11 + 55 = 66 DERIVED=np.stack([ # entry-0037-style derived set MIRf[ti,hi,wi], TIRf[ti,hi,wi], doz[ti,hi,wi], mir_z[ti,hi,wi], doz_z[ti,hi,wi], tdiff_z[ti,hi,wi], bg_mu[ti,hi,wi], bg_sd[ti,hi,wi], (MIRf[ti,hi,wi]-bg_mu[ti,hi,wi])/bg_sd[ti,hi,wi], ],1) COMBINED=np.concatenate([DERIVED,SPECTRAL],1) clat=LAT[ti,hi,wi]; clon=LON[ti,hi,wi] groups=(np.floor(clat/0.5).astype(int)*1000+np.floor(clon/0.5).astype(int)) # 0.5deg blocks y=labels.astype(int) def cv_auc(X,y,groups,shuffle=False,seed=0): yy=y.copy() if shuffle: rng=np.random.RandomState(seed); yy=rng.permutation(yy) uniq=np.unique(groups); k=min(5,len(uniq),int(yy.sum())) if k<2 or yy.sum()<5: return None,None,k gkf=GroupKFold(n_splits=k); aucs=[] for tr,te in gkf.split(X,yy,groups): if yy[tr].sum()<2 or yy[te].sum()<1 or (1-yy[te]).sum()<1: continue clf=HistGradientBoostingClassifier(max_iter=300,learning_rate=0.06,max_depth=4, l2_regularization=1.0,random_state=seed) clf.fit(X[tr],yy[tr]); p=clf.predict_proba(X[te])[:,1] try: aucs.append(roc_auc_score(yy[te],p)) except: pass if not aucs: return None,None,k return float(np.mean(aucs)),float(np.std(aucs)),k res={"n_frames":int(T),"date_range":[str(ftime[0]),str(ftime[-1])], "n_candidates":int(nC),"n_real":int(y.sum()),"n_false":int((1-y).sum()), "real_frac":float(y.mean()),"threshold":THRESH,"channels":CH} for name,X in [("DERIVED_0037baseline",DERIVED),("SPECTRAL_full_raw",SPECTRAL),("COMBINED",COMBINED)]: m,s,k=cv_auc(X,y,groups); res[name]={"cv_auc":m,"cv_std":s,"folds":k,"n_feat":X.shape[1]} log("%s spatial-CV AUC=%s +/- %s (k=%s, feat=%d)"%(name,m,s,k,X.shape[1])) # placebo on the strongest set pm,ps,pk=cv_auc(SPECTRAL,y,groups,shuffle=True,seed=7) res["PLACEBO_spectral_shuffled"]={"cv_auc":pm,"cv_std":ps,"folds":pk} log("PLACEBO (shuffled labels) AUC=%s (should be ~0.5)"%pm) json.dump(res,open(OUT,'w'),indent=1) np.savez_compressed(FEAT_OUT,derived=DERIVED,spectral=SPECTRAL,y=y,lat=clat,lon=clon, groups=groups,feat_names_spectral=np.array(CH+["%s-%s"%(CH[a],CH[b]) for (a,b) in pairs])) log("WROTE %s and %s"%(OUT,FEAT_OUT)) print("RESULT_JSON_BEGIN");print(json.dumps(res,indent=1));print("RESULT_JSON_END")