#!/usr/bin/env python3 """Accounting for the candidates. 1. Catalogue accounting of every candidate at f >= 0.1 (lines, UV, radio, AGN colour, post-starburst, photometry flags) and the mid-IR / radio test (q22) with FIRST and NVSS. 2. Legacy Surveys image cut-outs of the 'clean' candidates (no catalogue explanation) at f >= 0.1, in contact sheets (RAW/cutouts, figures/supp_clean_candidates.png). 3. Merge with code/vetting_manual.csv (visual classes entered by hand after looking at the contact sheets; kept in the repository so that the pipeline is re-runnable) and recompute the 'clean' limits with the visually explained objects removed -> data/limits_vetted.csv """ import json import os import urllib.parse import urllib.request import numpy as np import pandas as pd import mirq_common as M import s05_limits as L CUT = "https://www.legacysurvey.org/viewer/cutout.jpg" CUTDIR = os.path.join(M.RAW, "cutouts") MANUAL = os.path.join(M.HERE, "vetting_manual.csv") # Visual classes (entered by hand in vetting_manual.csv after looking at the contact sheets): # EXPLAINED positive image evidence for a conventional origin of the mid-IR flux or for a measurement failure # POSSIBLE a plausible but unproven conventional origin (disc morphology, very close companion, ...); # these are NOT removed from the 'unexplained' counts # everything else ('smooth_early_type') shows nothing unusual in the images EXPLAINED_CLASSES = {"artifact", "blend_neighbour", "merger_tidal", "spiral_disc"} POSSIBLE_CLASSES = {"blue_disc", "inclined_disc", "close_pair", "near_bright_star", "bad_image"} FIRST_LIM, NVSS_LIM = 1.0, 2.5 # mJy, catalogue detection limits N_SAMPLE_F01 = 48 # random clean candidates at (f=0.1, T=300 K) inspected for characterisation DATALAB = "https://datalab.noirlab.edu/tap/sync" IRSA = "https://irsa.ipac.caltech.edu/TAP/sync" _IRSA_UP = None def irsa_up(): global _IRSA_UP if _IRSA_UP is None: try: with urllib.request.urlopen(IRSA + "?" + urllib.parse.urlencode(dict( QUERY="SELECT TOP 1 ra FROM allwise_p3as_psd", FORMAT="CSV")), timeout=60) as r: _IRSA_UP = r.read().startswith(b"ra") except Exception: # noqa _IRSA_UP = False return _IRSA_UP def neighbours(vid, ra, dec, radius=25.0): """Mid-IR sources near a candidate: (a) AllWISE sources and (b) Legacy Surveys DR9 Tractor sources with their forced W3/W4 fluxes, within `radius` arcsec. Cached.""" os.makedirs(os.path.join(M.RAW, "neigh"), exist_ok=True) r = radius / 3600.0 qa = ("SELECT designation,ra,dec,w1mpro,w2mpro,w3mpro,w4mpro,w1snr,w3snr,w4snr,cc_flags,ext_flg FROM allwise_p3as_psd " "WHERE CONTAINS(POINT('ICRS',ra,dec),CIRCLE('ICRS',%.6f,%.6f,%.6f))=1" % (ra, dec, r)) def chk(b): if not b.startswith(b"designation"): raise RuntimeError("bad payload %r" % b[:100]) try: # IRSA is optional here (it was intermittently unavailable); the Legacy Surveys list is enough if not irsa_up(): raise RuntimeError("IRSA unavailable") fa = M.fetch(IRSA + "?" + urllib.parse.urlencode(dict(QUERY=qa, FORMAT="CSV")), os.path.join(M.RAW, "neigh", "%s_allwise.csv" % vid), min_bytes=20, timeout=120, tries=2, check=chk) except Exception: # noqa fa = None ql = ("SELECT ra,dec,type,flux_r,flux_w1,flux_w3,flux_w4,flux_ivar_w3,flux_ivar_w4,shape_r FROM ls_dr9.tractor " "WHERE 't'=q3c_radial_query(ra,dec,%.6f,%.6f,%.6f)" % (ra, dec, r)) def chk2(b): if not b.startswith(b"ra,dec"): raise RuntimeError("bad payload %r" % b[:100]) fl_ = M.fetch(DATALAB, os.path.join(M.RAW, "neigh", "%s_ls.csv" % vid), min_bytes=10, timeout=300, check=chk2, data=urllib.parse.urlencode(dict(REQUEST="doQuery", LANG="ADQL", FORMAT="csv", QUERY=ql)).encode()) l = pd.read_csv(fl_) l["sep"] = M.angsep_arcsec(ra, dec, l.ra.values, l["dec"].values) out = dict(n_ls=int(len(l))) if fa is not None: a = pd.read_csv(fa) a["sep"] = M.angsep_arcsec(ra, dec, a.ra.values, a["dec"].values) oth = a[a.sep > 3.0] out["n_allwise"] = int(len(a)) out["allwise_other_n"] = int(len(oth)) out["allwise_other_max_w4snr"] = float(oth.w4snr.max()) if len(oth) else np.nan out["allwise_other_sep_of_max_w4"] = float(oth.sep[oth.w4snr.idxmax()]) if len(oth) and oth.w4snr.notna().any() else np.nan lo = l[(l.sep > 1.5) & (l.type != "DUP")] if len(lo): k = lo.flux_w4.idxmax() out["ls_other_max_w4_mJy"] = float(M.ls_nmgy_to_wise_jy(lo.flux_w4[k], 4) * 1e3) out["ls_other_sep_of_max_w4"] = float(lo.sep[k]) k3 = lo.flux_w3.idxmax() out["ls_other_max_w3_mJy"] = float(M.ls_nmgy_to_wise_jy(lo.flux_w3[k3], 3) * 1e3) out["ls_other_sep_of_max_w3"] = float(lo.sep[k3]) return out def get_cutout(ra, dec, name, pixscale=0.4, size=200, layers=("ls-dr9", "ls-dr10", "sdss"), suffix=""): """Cached JPEG cut-out from the Legacy Surveys viewer; the first layer that returns an image is used.""" os.makedirs(CUTDIR, exist_ok=True) for layer in layers: dest = os.path.join(CUTDIR, "%s%s_%s.jpg" % (name, suffix, layer)) if not (os.path.exists(dest) and os.path.getsize(dest) > 2500): url = CUT + "?" + urllib.parse.urlencode(dict(ra="%.5f" % ra, dec="%.5f" % dec, layer=layer, pixscale=pixscale, size=size)) try: with urllib.request.urlopen(url, timeout=90) as r: blob = r.read() with open(dest, "wb") as fh: fh.write(blob) except Exception: # noqa continue if os.path.getsize(dest) > 2500: return dest return None # the three views of each inspected object: (suffix, pixscale ["/px], size [px], layers, panel label) VIEWS = (("", 0.4, 200, ("ls-dr9", "ls-dr10", "sdss"), "optical 80\""), ("_zoom", 0.15, 200, ("ls-dr9", "ls-dr10", "sdss"), "optical 30\""), ("_w34", 1.375, 128, ("unwise-w3w4",), "unWISE W3/W4 176\"")) def contact_sheets(c, tag, nobj_row=3, rows_per_sheet=6, out_fig=None): """Three panels per object: optical 80", optical zoom 30", unWISE W3/W4 176" (north up, east left).""" import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import matplotlib.image as mpimg paths = [] per_sheet = nobj_row * rows_per_sheet for k0 in range(0, len(c), per_sheet): sub = c.iloc[k0:k0 + per_sheet] nrow = int(np.ceil(len(sub) / nobj_row)) fig, axes = plt.subplots(nrow, 3 * nobj_row, figsize=(2.0 * 3 * nobj_row, 2.3 * nrow), squeeze=False) for ax in axes.ravel(): ax.axis("off") for k, (_, r) in enumerate(sub.iterrows()): for j, (suf, ps, size, layers, lab) in enumerate(VIEWS): ax = axes[k // nobj_row, 3 * (k % nobj_row) + j] fn = get_cutout(r.ra, r["dec"], "g%d" % r.id, pixscale=ps, size=size, layers=layers, suffix=suf) if fn: try: im = mpimg.imread(fn) ax.imshow(im) if suf == "_w34": # mark the galaxy position (the WISE W4 beam is 12" FWHM) ax.plot([im.shape[1] / 2.0], [im.shape[0] / 2.0], marker="o", ms=9, mfc="none", mec="#ffffff", mew=0.7) except Exception: # noqa pass if j == 0: ax.set_title("%s z=%.3f f>=%.1f %s\n%s" % (r.vid, r.z, r.f_max, r.line_class, r.clean_ft_list[:38]), fontsize=6.3, loc="left") elif j == 1: ax.set_title("W1-W2=%.2f W1-W3=%.1f\nF3=%.2f+-%.2f mJy" % (r.w12, r.w13, r.F3, r.e3), fontsize=6.3, loc="left") else: ax.set_title("F4=%.1f+-%.1f mJy (unWISE W3/W4)\n%s" % (r.F4, r.e4, r.visual_class), fontsize=6.3, loc="left") fig.tight_layout(pad=0.3) p = os.path.join(CUTDIR, "sheet_%s_%02d.png" % (tag, k0 // per_sheet)) fig.savefig(p, dpi=110) if out_fig and k0 == 0: fig.savefig(out_fig, dpi=110) plt.close(fig) paths.append(p) return paths def main(make_sheets=True, fetch_neighbours=True): d = L.load_master() fl = L.flags(d) eng = L.Engine(d) summ = {} # ---------- every candidate at f >= 0.1 (any T of the grid), all populations ---------- rows = {} for T in M.T_GRID: for f in (0.1, 0.3): use_col = eng.predicted_w12(f, T) < 0.3 clean = L.clean_mask(fl, use_agncol=use_col) _, cand, _ = eng.run(f, T) for i in np.where(cand)[0]: e = rows.setdefault(i, dict(f_max=0.0, ft=[], clean_ft=[])) e["f_max"] = max(e["f_max"], f) e["ft"].append("%.1f@%d" % (f, T)) if clean[i]: e["clean_ft"].append("%.1f@%d" % (f, T)) idx = np.array(sorted(rows)) c = d.loc[idx].copy() c["f_max"] = [rows[i]["f_max"] for i in idx] c["ft_list"] = [";".join(rows[i]["ft"]) for i in idx] c["clean_ft_list"] = [";".join(rows[i]["clean_ft"]) for i in idx] c["is_clean_cand"] = c.clean_ft_list != "" for k in fl.columns: c[k] = fl.loc[idx, k].values c["w12"] = c.w1mpro - c.w2mpro with np.errstate(invalid="ignore", divide="ignore"): c["w13"] = 2.5 * np.log10(np.clip(c.F3 / c.F1, 1e-3, None) * M.WISE_F0[1] / M.WISE_F0[3]) c["w14"] = 2.5 * np.log10(np.clip(c.F4 / c.F1, 1e-3, None) * M.WISE_F0[1] / M.WISE_F0[4]) # mid-IR / radio: q22 = log10(S_22um / S_1.4GHz); FIRST preferred, NVSS otherwise, lower limit if neither s14 = c.first_fint10.where(c.first_fint10.notna(), c.nvss_s14) c["s14_mJy"] = s14 c["q22"] = np.log10(c.F4 / s14) c["q22_lowlim"] = np.where(s14.isna() & (c.F4 > 3 * c.e4), np.log10(np.clip(c.F4, 1e-3, None) / FIRST_LIM), np.nan) c = c.sort_values(["is_clean_cand", "f_max", "w13"], ascending=[False, False, False]).reset_index(drop=True) c["vid"] = ["C%03d" % k for k in range(len(c))] summ["n_cand_f>=0.1"] = int(len(c)) summ["n_cand_f>=0.3"] = int((c.f_max >= 0.3).sum()) summ["n_clean_cand_f>=0.1"] = int(c.is_clean_cand.sum()) for k in fl.columns: summ["cand_f>=0.1_" + k] = int(c[k].sum()) summ["cand_f>=0.1_no_flag_at_all"] = int((~c[list(fl.columns)].any(axis=1)).sum()) det = c.s14_mJy.notna() & (c.F4 > 3 * c.e4) summ["radio"] = dict(n_with_radio=int(c.s14_mJy.notna().sum()), n_q22_measured=int(det.sum()), q22_median=float(c.q22[det].median()) if det.sum() else None, n_q22_gt_1p5=int((c.q22[det] > 1.5).sum()), n_q22_gt_2=int((c.q22[det] > 2.0).sum()), n_q22_lt_0=int((c.q22[det] < 0).sum()), n_lowlim=int(c.q22_lowlim.notna().sum()), n_lowlim_gt_1p5=int((c.q22_lowlim > 1.5).sum()), n_lowlim_gt_2=int((c.q22_lowlim > 2.0).sum()), n_lowlim_gt_2_clean=int(((c.q22_lowlim > 2.0) & c.is_clean_cand).sum()), lowlim_median=float(c.q22_lowlim.median()) if c.q22_lowlim.notna().any() else None) # ---------- manual visual classes ---------- if os.path.exists(MANUAL): man = pd.read_csv(MANUAL, dtype={"id": "int64"}) c = c.merge(man[["id", "visual_class", "note"]], on="id", how="left") else: c["visual_class"], c["note"] = np.nan, np.nan c["visual_class"] = c.visual_class.fillna("not_inspected") c.loc[~c.is_clean_cand & (c.visual_class == "not_inspected"), "visual_class"] = "n/a (catalogue-explained)" keep = ["vid", "id", "ra", "dec", "z", "f_max", "ft_list", "clean_ft_list", "is_clean_cand", "w12", "w13", "w14", "F1", "F3", "F4", "e3", "e4", "L1", "L3", "L4", "ext_flg", "theta_e", "line_class", "ew_ha", "d4000", "lick_hd_a", "nuv_r", "s14_mJy", "q22", "q22_lowlim", "veldisp", "Mr", "x_line", "x_agncol", "x_uv", "x_psb", "x_radio", "x_phot", "visual_class", "note"] for k in ("ra", "dec"): # keep full positional precision in the CSV outputs c[k + "_deg"] = c[k].map(lambda v: "%.5f" % v) keep = [k + "_deg" if k in ("ra", "dec") else k for k in keep] c[keep].to_csv(os.path.join(M.DATA, "candidates_f0.1plus.csv"), index=False, float_format="%.5g") cc = c[c.is_clean_cand].reset_index(drop=True) # inspection sets: A = every clean candidate at f = 0.3 (any T); B = seeded random sample of the clean # candidates at (f = 0.1, T = 300 K), for characterisation only cc["set_A"] = cc.clean_ft_list.str.contains("0.3@") b_pool = cc.index[cc.clean_ft_list.str.contains("0.1@300") & ~cc.set_A].values # Set B was drawn once (seeded, 48 of the pool) in the first pass and inspected by hand; the draw is frozen # in vetting_manual.csv (inspection_set == 'B') so that later changes of the clean sample (revision: failed # H-alpha fits are no longer 'lineless') do not silently replace inspected objects by uninspected ones. man_b = set() if os.path.exists(MANUAL): mm = pd.read_csv(MANUAL, dtype={"id": "int64"}) man_b = set(mm.id[mm.inspection_set == "B"].values) if man_b: b_pick = np.array([i for i in b_pool if cc.id[i] in man_b]) summ["n_setB_frozen_draw"] = int(len(man_b)) else: rng = np.random.default_rng(L.RNG_SEED) b_pick = np.sort(rng.choice(b_pool, size=min(N_SAMPLE_F01, len(b_pool)), replace=False)) cc["set_B"] = False cc.loc[b_pick, "set_B"] = True summ["n_clean_cand_f0.3_setA"] = int(cc.set_A.sum()) summ["n_clean_cand_f0.1_T300_pool"] = int(len(b_pool)) summ["n_setB_sampled"] = int(cc.set_B.sum()) # neighbours of the set-A objects nb_rows = [] for _, r in (cc[cc.set_A].iterrows() if fetch_neighbours else ()): try: nb = neighbours("g%d" % r.id, r.ra, r["dec"]) except Exception as e: # noqa nb = dict(error=str(e)[:80]) nb["id"] = r.id nb_rows.append(nb) if nb_rows: cc = cc.merge(pd.DataFrame(nb_rows), on="id", how="left") def tier(v): return ("explained" if v in EXPLAINED_CLASSES else "possible" if v in POSSIBLE_CLASSES else "not_inspected" if v == "not_inspected" else "nothing_unusual") cc["tier"] = cc.visual_class.map(tier) for nm in ("A", "B"): sub = cc[cc["set_" + nm]] summ["set%s_visual_class_counts" % nm] = sub.visual_class.value_counts().to_dict() summ["set%s_tier_counts" % nm] = sub.tier.value_counts().to_dict() extra = [k for k in cc.columns if k.startswith(("allwise_other", "ls_other", "n_allwise", "n_ls"))] cc[keep + ["set_A", "set_B", "tier"] + extra].to_csv(os.path.join(M.DATA, "candidates_clean.csv"), index=False, float_format="%.5g") # catalogue properties of the clean candidates at (f = 0.1, T = 300 K) against the whole clean sample clean_all = d[L.clean_mask(fl, True) & eng.P["ok"]] pool = cc[cc.clean_ft_list.str.contains("0.1@300")] summ["clean_cand_f0.1_T300_vs_clean_sample_medians"] = { k: [float(pool[k].median()), float(clean_all[k].median())] for k in ("z", "ba", "gr", "d4000", "veldisp", "sersic", "ew_ha", "Mr")} summ["clean_cand_f0.1_T300_frac_retired"] = [float((pool.line_class == "retired").mean()), float((clean_all.line_class == "retired").mean())] summ["clean_cand_f0.1_T300_frac_ba_lt_0.5"] = [float((pool.ba < 0.5).mean()), float((clean_all.ba < 0.5).mean())] if make_sheets and len(cc): contact_sheets(cc[cc.set_A], "setA", out_fig=os.path.join(M.FIGS, "supp_clean_candidates.png")) contact_sheets(cc[cc.set_B], "setB") # ---------- limits after visual vetting of the clean candidates ---------- vis_expl = set(cc.id[cc.visual_class.isin(EXPLAINED_CLASSES)].values) vis_poss = set(cc.id[cc.visual_class.isin(POSSIBLE_CLASSES)].values) is_insp = d.id.isin(set(cc.id[cc.visual_class != "not_inspected"].values)).values is_vis = d.id.isin(vis_expl).values is_poss = d.id.isin(vis_poss).values rows = [] for T in M.T_GRID: for f in M.F_VAR: use_col = eng.predicted_w12(f, T) < 0.3 clean = L.clean_mask(fl, use_agncol=use_col) r, cand, test = eng.run(f, T, clean) n_un = int((cand & ~is_vis).sum()) n_ni = int((cand & ~is_insp).sum()) n_po = int((cand & is_poss).sum()) # 'nothing unusual' count is meaningful only where every candidate was inspected n_nu = (n_un - n_po) if n_ni == 0 else np.nan rows.append(dict(T=T, f=f, n_test_clean=r["n_test"], n_cand_clean=r["n_cand"], n_cand_clean_visually_explained=r["n_cand"] - n_un, n_cand_unexplained=n_un, n_cand_possible=n_po, n_cand_nothing_unusual=n_nu, n_cand_not_inspected=n_ni, ul_clean=L.limit(r["n_cand"], r["n_test"], r["recov_eff"]), ul_unexplained=L.limit(n_un, r["n_test"], r["recov_eff"]), ul_nothing_unusual=(L.limit(int(n_nu), r["n_test"], r["recov_eff"]) if n_ni == 0 else np.nan), ul_zero_clean=L.limit(0, r["n_test"], r["recov_eff"]))) v = pd.DataFrame(rows) v.to_csv(os.path.join(M.DATA, "limits_vetted.csv"), index=False, float_format="%.5g") print(v.to_string(index=False, float_format=lambda x: "%.3g" % x)) # the mid-IR-bright, radio-undetected objects (q22 lower limit > 2, the outlier definition of Garrett 2015) hi = c[c.q22_lowlim > 2.0] summ["q22_lowlim_gt_2_objects"] = [dict(vid=r.vid, id=int(r.id), F4_mJy=float(r.F4), F3_mJy=float(r.F3), cc_flags=str(r.cc_flags), line_class=r.line_class, w12=float(r.w12), is_clean_cand=bool(r.is_clean_cand), visual_class=r.visual_class) for _, r in hi.iterrows()] json.dump(summ, open(os.path.join(M.DATA, "vetting_summary.json"), "w"), indent=1) print(json.dumps(summ, indent=1)) if __name__ == "__main__": main()