"""Step 7: image-cutout check of every excess star (AllWISE Atlas images through CDS hips2fits). A swarm (or a disc) radiates at the position of the star. A background galaxy that fell into the 6.5 arcsec (W3) / 12 arcsec (W4) beam does not: it pulls the W3/W4 light centroid away from the position of the star. For every excess star and for a comparison ("null") sample of non-excess stars matched three-to-one in catalogue W3 (W4) signal-to-noise we * download 88 x 88 arcsec cutouts in W1, W3, W4 centred on the Gaia position of the star propagated to the WISE epoch (and, for excess stars, an optical r-band cutout: Pan-STARRS1 north of Dec -30, DESI Legacy DR10 or DSS2-red elsewhere); * measure the Gaussian-windowed centroid in W3 and W4 and its offset d from the star; * calibrate the offset noise sigma_d = sqrt(s0^2 + (k / SNR_cat)^2) on the null stars and flag excess stars whose normalised offset d / sigma_d exceeds the 99th percentile of the null stars (so 1% of genuine on-star sources are flagged by construction), or whose centroid cannot be measured at all because it walks off to a neighbouring source more than 12 arcsec away (0.2% of the null stars). Outputs data/image_centroids.csv : per excess star, offsets, normalised offsets, flag data/image_centroid_null.csv : the same for the comparison stars data/image_centroid_calib.json : s0, k per band and the comparison-sample exceedance rate RAW/cutouts/*.fits, RAW/sheets/*.png (contact sheets for visual inspection) """ import json import os import sys import urllib.parse from concurrent.futures import ThreadPoolExecutor import numpy as np import pandas as pd from astropy.io import fits from common import RAW, DATA, http_get import importlib inj = importlib.import_module("04_injection") vq = importlib.import_module("06_vet_queries") CUT = os.path.join(RAW, "cutouts") SHEETS = os.path.join(RAW, "sheets") HIPS2FITS = "https://alasky.cds.unistra.fr/hips-image-services/hips2fits" NPIX, PIXSCALE = 64, 1.375 # WISE cutouts: 88 arcsec, native Atlas pixel scale OPT_NPIX, OPT_PIX = 120, 0.5 # optical cutouts: 60 arcsec HIPS = {"W1": "CDS/P/allWISE/W1", "W3": "CDS/P/allWISE/W3", "W4": "CDS/P/allWISE/W4", "ps1r": "CDS/P/PanSTARRS/DR1/r", "lsr": "CDS/P/DESI-Legacy-Surveys/DR10/r", "dss": "CDS/P/DSS2/red"} SIG_WIN = {"W1": 3.6, "W3": 3.6, "W4": 7.0} # arcsec, Gaussian window sigma (about the Atlas PSF sigma) N_NULL = 3 # null stars per excess star and band NULL_PCT = 99.0 # percentile of the null distribution used as the flag threshold MIN_SNR_IMG = 4.0 # window S/N needed for a centroid to be used def fetch(sid, key, ra, dec): os.makedirs(CUT, exist_ok=True) out = os.path.join(CUT, f"{sid}_{key}.fits") if os.path.exists(out) and os.path.getsize(out) > 2880: return out wise = key in ("W1", "W3", "W4") n, pix = (NPIX, PIXSCALE) if wise else (OPT_NPIX, OPT_PIX) q = urllib.parse.urlencode(dict(hips=HIPS[key], width=n, height=n, fov=n * pix / 3600.0, projection="TAN", coordsys="icrs", ra=f"{ra:.6f}", dec=f"{dec:.6f}", format="fits")) try: data = http_get(HIPS2FITS + "?" + q, timeout=120) except Exception as exc: # noqa: BLE001 (a missing cutout is recorded as missing, not fatal) print(f"cutout failed {sid} {key}: {exc}", flush=True) return None if isinstance(data, str): data = data.encode("latin-1") with open(out, "wb") as f: f.write(data) return out def read_img(sid, key, with_centre=False): """Image array (and the 0-based pixel coordinates of the requested position).""" p = os.path.join(CUT, f"{sid}_{key}.fits") if not os.path.exists(p): return (None, None) if with_centre else None try: with fits.open(p) as h: img = np.asarray(h[0].data, dtype=float) c0 = (float(h[0].header["CRPIX1"]) - 1.0, float(h[0].header["CRPIX2"]) - 1.0) except Exception: # noqa: BLE001 return (None, None) if with_centre else None if img.ndim != 2 or not np.isfinite(img).any() or np.nanstd(img) == 0: return (None, None) if with_centre else None return (img, c0) if with_centre else img def windowed_centroid(img, c0, sig_arcsec, pix=PIXSCALE): """Gaussian-windowed centroid started at the star's position c0 = (x0, y0) [pixels]. Returns (dx, dy) in arcsec relative to the star (dx positive = +x pixel direction = west) and the window S/N; NaN if the centroid wanders off (> 12 arcsec). """ n = img.shape[0] yy, xx = np.mgrid[0:n, 0:n] cx, cy = c0 r = np.hypot(xx - cx, yy - cy) * pix ann = (r > 25) & (r < 42) & np.isfinite(img) if ann.sum() < 50: return np.nan, np.nan, np.nan bkg = np.median(img[ann]) noise = 1.4826 * np.median(np.abs(img[ann] - bkg)) im = np.where(np.isfinite(img), img - bkg, 0.0) sw = sig_arcsec / pix x, y = cx, cy for _ in range(30): w = np.exp(-0.5 * ((xx - x) ** 2 + (yy - y) ** 2) / sw ** 2) w[np.hypot(xx - x, yy - y) > 3 * sw] = 0 tot = (w * im).sum() if not np.isfinite(tot) or tot <= 0: return np.nan, np.nan, 0.0 xn = x + 2 * ((w * im * xx).sum() / tot - x) yn = y + 2 * ((w * im * yy).sum() / tot - y) if np.hypot(xn - cx, yn - cy) * pix > 12: return np.nan, np.nan, tot / (noise * np.sqrt((w ** 2).sum()) + 1e-30) done = np.hypot(xn - x, yn - y) < 0.01 x, y = xn, yn if done: break w = np.exp(-0.5 * ((xx - x) ** 2 + (yy - y) ** 2) / sw ** 2) snr = (w * im).sum() / (noise * np.sqrt((w ** 2).sum()) + 1e-30) return (x - cx) * pix, (y - cy) * pix, snr def comparison_sample(d, n_null=N_NULL): """n_null non-excess stars per excess star and band, nearest in catalogue W3 (W4) S/N.""" rows = [] used = set() for b in ("3", "4"): det = d[d["det" + b]] pool = d[d["q" + b] & ~d.det & (d["chi" + b].abs() < 2) & (d[f"w{b}snr"] > 3)].sort_values(f"w{b}snr") ps = np.log10(pool[f"w{b}snr"].values) psid = pool.sid.values for _, s in det.sort_values("sid").iterrows(): target = np.log10(max(s[f"w{b}snr"], 3.0)) j = int(np.clip(np.searchsorted(ps, target), 0, len(ps) - 1)) got = 0 for step in range(0, 5000): # nearest unused for jj in (j + step, j - step): if got < n_null and 0 <= jj < len(ps) and psid[jj] not in used: used.add(psid[jj]) rows.append(dict(sid=psid[jj], band="W" + b, match_sid=s.sid)) got += 1 if got >= n_null: break return pd.DataFrame(rows) def download(d, comp): det = d[d.det] pos = vq.positions(det).set_index("sid") jobs = [] for sid, p in pos.iterrows(): for k in ("W1", "W3", "W4"): jobs.append((sid, k, p.ra_wise, p.dec_wise)) jobs.append((sid, "ps1r" if p.dec_wise > -29.5 else "lsr", p.ra_wise, p.dec_wise)) cst = d[d.sid.isin(comp.sid)] cpos = vq.positions(cst).set_index("sid") for _, r in comp.iterrows(): p = cpos.loc[r.sid] jobs.append((r.sid, r.band, p.ra_wise, p.dec_wise)) jobs = [j for j in jobs if not os.path.exists(os.path.join(CUT, f"{j[0]}_{j[1]}.fits"))] print(f"downloading {len(jobs)} cutouts", flush=True) with ThreadPoolExecutor(4) as ex: # network I/O only list(ex.map(lambda j: fetch(*j), jobs)) # southern stars outside the Legacy footprint: fall back to DSS2 red fb = [(sid, "dss", p.ra_wise, p.dec_wise) for sid, p in pos.iterrows() if p.dec_wise <= -29.5 and read_img(sid, "lsr") is None and not os.path.exists(os.path.join(CUT, f"{sid}_dss.fits"))] with ThreadPoolExecutor(4) as ex: list(ex.map(lambda j: fetch(*j), fb)) def measure(sids, bands): rows = [] for sid in sids: r = dict(sid=sid) for b in bands: img, c0 = read_img(sid, b, with_centre=True) if img is None: r.update({f"dx{b}": np.nan, f"dy{b}": np.nan, f"snr_img{b}": np.nan}) continue dx, dy, snr = windowed_centroid(img, c0, SIG_WIN[b]) r.update({f"dx{b}": dx, f"dy{b}": dy, f"snr_img{b}": snr}) rows.append(r) return pd.DataFrame(rows) def calibrate_null(null, b): """Zero point (median dx, dy), noise model sigma_d = sqrt(s0^2 + (k/SNR_cat)^2) and the empirical NULL_PCT percentile of d / sigma_d, all from the non-excess comparison stars.""" v = null[np.isfinite(null["dx" + b]) & (null["snr_img" + b] > MIN_SNR_IMG)] x0, y0 = float(np.median(v["dx" + b])), float(np.median(v["dy" + b])) off = np.hypot(v["dx" + b] - x0, v["dy" + b] - y0).values snr = v["snr_cat"].values med_rayleigh = np.sqrt(2 * np.log(2)) hi = snr > np.quantile(snr, 0.75) s0 = np.median(off[hi]) / med_rayleigh lo = snr < np.quantile(snr, 0.5) ks = np.linspace(0.0, 30, 601) m = np.array([np.median(off[lo] / np.sqrt(s0 ** 2 + (k / snr[lo]) ** 2)) for k in ks]) k = float(ks[np.argmin(np.abs(m - med_rayleigh))]) dn = off / np.sqrt(s0 ** 2 + (k / snr) ** 2) return dict(x0=x0, y0=y0, s0=float(s0), k=k, n_null=int(len(null)), n_null_usable=int(len(v)), median_norm=float(np.median(dn)), threshold=float(np.percentile(dn, NULL_PCT)), rayleigh_99=3.035, frac_null_above_rayleigh99=float((dn > 3.035).mean())) def normalise(t, b, cal): off = np.hypot(t["dx" + b] - cal["x0"], t["dy" + b] - cal["y0"]) t["off" + b] = off t["n" + b] = off / np.sqrt(cal["s0"] ** 2 + (cal["k"] / t["snr_cat" + ("" if "snr_cat" in t else b)].clip(lower=1.0)) ** 2) t.loc[~(t["snr_img" + b] > MIN_SNR_IMG), "n" + b] = np.nan return t def sheets(d, cen, order, tag, per=8, outdir=None, labels=None, dpi=110): """Contact sheets: optical r | W1 | W3 | W4, 60 arcsec on a side, north up, east left. Red cross and 6 arcsec circle: the star at the WISE epoch; blue x: W3/W4 light centroid.""" import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt outdir = outdir or SHEETS os.makedirs(outdir, exist_ok=True) det = d.set_index("sid") cen = cen.set_index("sid") out = [] for k0 in range(0, len(order), per): chunk = order[k0:k0 + per] fig, axs = plt.subplots(len(chunk), 4, figsize=(8.4, 2.1 * len(chunk)), squeeze=False) for i, sid in enumerate(chunk): s = det.loc[sid] c = cen.loc[sid] for j, key in enumerate(("opt", "W1", "W3", "W4")): ax = axs[i, j] ax.set_xticks([]) ax.set_yticks([]) if key == "opt": img, lab, c0 = None, "", None for kk in ("ps1r", "lsr", "dss"): img, c0 = read_img(sid, kk, with_centre=True) if img is not None: lab = kk break pix = OPT_PIX else: (img, c0), lab, pix = read_img(sid, key, with_centre=True), key, PIXSCALE if img is None: ax.text(0.5, 0.5, "no image", ha="center", va="center", transform=ax.transAxes) continue n = img.shape[0] fin = img[np.isfinite(img)] lo, hi = np.percentile(fin, [20, 99.7 if key == "opt" else 99.9]) # east to the left: x pixel index increases to the west ext = [(c0[0] + 0.5) * pix, (c0[0] + 0.5 - n) * pix, (-c0[1] - 0.5) * pix, (n - c0[1] - 0.5) * pix] ax.imshow(np.arcsinh((np.nan_to_num(img, nan=lo) - lo) / max(hi - lo, 1e-30) * 10), origin="lower", cmap="gray_r", extent=ext) ax.set_xlim(30, -30) ax.set_ylim(-30, 30) ax.plot([0], [0], "+", color="red", ms=14, mew=0.8) th = np.linspace(0, 2 * np.pi, 100) ax.plot(6 * np.cos(th), 6 * np.sin(th), color="red", lw=0.5) if key in ("W3", "W4") and np.isfinite(c.get("dx" + key, np.nan)): # image x increases to the west (RA decreases): dx (pixels, +x) -> east offset = -dx ax.plot([-c["dx" + key]], [c["dy" + key]], "x", color="blue", ms=7, mew=0.8) ax.set_title(lab, fontsize=7, pad=1) if j == 0: head = labels[sid] if labels and sid in labels else str(sid)[-7:] ax.set_ylabel(f"{head}\n{s.cls} d={s.dist_pc:.0f}pc\nG={s.Gmag:.1f}", fontsize=6) if j == 2: ax.set_xlabel(f"chi3={s.chi3:.1f} off={c.get('offW3', np.nan):.1f}\" n={c.get('nW3', np.nan):.1f}", fontsize=6, labelpad=1) if j == 3: ax.set_xlabel(f"chi4={s.chi4:.1f} off={c.get('offW4', np.nan):.1f}\" n={c.get('nW4', np.nan):.1f}", fontsize=6, labelpad=1) if j == 1: ax.set_xlabel(str(s.get("explanation", ""))[:28], fontsize=6, labelpad=1) fig.subplots_adjust(left=0.08, right=0.995, top=0.985, bottom=0.02, wspace=0.03, hspace=0.22) p = os.path.join(outdir, f"{tag}_{k0 // per:02d}.png" if len(order) > per else f"{tag}.png") fig.savefig(p, dpi=dpi) plt.close(fig) out.append(p) return out def main(make_sheets=True): d = inj.load_stars() comp = comparison_sample(d) download(d, comp) det = d[d.det] cen = measure(det.sid.values, ("W1", "W3", "W4")) cen = cen.merge(det[["sid", "w3snr", "w4snr", "det3", "det4"]], on="sid") null, calib = {}, {} for b in ("W3", "W4"): nb = measure(comp[comp.band == b].sid.values, (b,)) nb["band"] = b nb = nb.merge(d[["sid", f"w{b[1]}snr"]].rename(columns={f"w{b[1]}snr": "snr_cat"}), on="sid") calib[b] = calibrate_null(nb, b) null[b] = normalise(nb, b, calib[b]) cen["snr_cat" + b] = cen[f"w{b[1]}snr"] cen = normalise(cen, b, calib[b]) # 'lost': the windowed centroid walked out of the 12 arcsec search region (to a brighter # neighbouring source) or found no peak at the star; this happens for ~0.2% of the null stars cen["lost_W3"] = cen.det3 & ~np.isfinite(cen.nW3) cen["lost_W4"] = cen.det4 & ~np.isfinite(cen.nW4) cen["flag_off_W3"] = ((cen.nW3 > calib["W3"]["threshold"]) & cen.det3) | cen.lost_W3 cen["flag_off_W4"] = ((cen.nW4 > calib["W4"]["threshold"]) & cen.det4) | cen.lost_W4 cen["flag_offset"] = cen.flag_off_W3 | cen.flag_off_W4 cen.drop(columns=["w3snr", "w4snr"]).to_csv(os.path.join(DATA, "image_centroids.csv"), index=False, float_format="%.4g") pd.concat(null.values()).to_csv(os.path.join(DATA, "image_centroid_null.csv"), index=False, float_format="%.4g") with open(os.path.join(DATA, "image_centroid_calib.json"), "w") as f: json.dump(calib, f, indent=1) print("centroid-noise calibration:", json.dumps(calib)) m = cen.merge(det[["sid", "cls"]], on="sid") for b, k in (("W3", "det3"), ("W4", "det4")): print(f"excess stars with a usable {b} centroid:", int((m[k] & np.isfinite(m["n" + b])).sum()), "of", int(m[k].sum()), "; offset-flagged:", int(m["flag_off_" + b].sum())) print(pd.crosstab(m.cls, m.flag_offset).to_string()) if make_sheets: cfile = os.path.join(DATA, "candidates.csv") if os.path.exists(cfile): c = pd.read_csv(cfile) dets = det.merge(c[["sid", "explanation"]], on="sid", how="left") order = c.sort_values(["explanation", "cls", "chi3"], ascending=[True, True, False]).sid.tolist() else: dets, order = det, det.sort_values("chi3", ascending=False).sid.tolist() sheets(dets, cen, order, "sheet") return cen if __name__ == "__main__": main(make_sheets="--no-sheets" not in sys.argv)