#!/usr/bin/env python3 """Vetting of the most under-luminous FP outliers. 1. AllWISE photometry for the complete SDSS and 6dF samples (CDS XMatch, 5 arcsec) -> mid-IR colours, and an r-W1 (J-W1) colour residual that flags optical-photometry problems / dust. 2. Image cut-outs (Legacy Surveys viewer) for every galaxy with dI <= CAND_THR, assembled in contact sheets. 3. Merges automatic diagnostics with the stored visual classification (code/vetting_manual.csv, entered by hand after looking at the contact sheets; kept under version control so the pipeline is re-runnable). """ import io import os import sys import time import urllib.parse import urllib.request import uuid import numpy as np import pandas as pd import fpk_common as F CAND_THR = -0.45 # dex: equivalent to >= 65 per cent of the light missing XMATCH = "http://cdsxmatch.u-strasbg.fr/xmatch/api/v1/sync" CUT = "https://www.legacysurvey.org/viewer/cutout.jpg" CUTDIR = os.path.join(F.RAW, "cutouts") def _multipart(fields, files): boundary = uuid.uuid4().hex body = io.BytesIO() for k, v in fields.items(): body.write(("--%s\r\nContent-Disposition: form-data; name=\"%s\"\r\n\r\n%s\r\n" % (boundary, k, v)).encode()) for k, (fn, blob) in files.items(): body.write(("--%s\r\nContent-Disposition: form-data; name=\"%s\"; filename=\"%s\"\r\n" "Content-Type: text/csv\r\n\r\n" % (boundary, k, fn)).encode()) body.write(blob) body.write(b"\r\n") body.write(("--%s--\r\n" % boundary).encode()) return body.getvalue(), "multipart/form-data; boundary=%s" % boundary def xmatch(df, cat2, dest, radius=5.0, chunk=40000): """Positional cross-match of df[id, ra, dec] with a VizieR catalogue through the CDS XMatch service.""" if os.path.exists(dest) and os.path.getsize(dest) > 200: return pd.read_csv(dest) parts = [] for i0 in range(0, len(df), chunk): sub = df.iloc[i0:i0 + chunk][["id", "ra", "dec"]] blob = sub.to_csv(index=False).encode() data, ctype = _multipart(dict(request="xmatch", distMaxArcsec=str(radius), RESPONSEFORMAT="csv", colRA1="ra", colDec1="dec", cat2=cat2), {"cat1": ("cat1.csv", blob)}) last = None for k in range(5): try: req = urllib.request.Request(XMATCH, data=data, headers={"Content-Type": ctype}) with urllib.request.urlopen(req, timeout=900) as r: txt = r.read().decode() parts.append(pd.read_csv(io.StringIO(txt))) break except Exception as e: # noqa last = e time.sleep(10 * (k + 1)) else: raise RuntimeError("xmatch failed: %s" % last) out = pd.concat(parts, ignore_index=True) out.to_csv(dest, index=False) return out def get_cutout(ra, dec, name, pixscale, layers=("ls-dr10", "ls-dr9", "sdss")): os.makedirs(CUTDIR, exist_ok=True) for layer in layers: dest = os.path.join(CUTDIR, "%s_%s.jpg" % (name, layer)) if not (os.path.exists(dest) and os.path.getsize(dest) > 1500): url = CUT + "?" + urllib.parse.urlencode(dict(ra="%.5f" % ra, dec="%.5f" % dec, layer=layer, pixscale=pixscale, size=200)) try: with urllib.request.urlopen(url, timeout=90) as r: blob = r.read() with open(dest, "wb") as fh: fh.write(blob) except Exception: continue try: from PIL import Image im = np.asarray(Image.open(dest).convert("RGB")) if im.std() > 4: # blank tiles (outside footprint) are uniform return im, layer except Exception: continue return None, None def contact_sheet(cands, fn, title, pixscale_col="pixscale"): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt n = len(cands) ncol = 6 nrow = int(np.ceil(n / ncol)) fig, axes = plt.subplots(nrow, ncol, figsize=(ncol * 2.3, nrow * 2.55)) axes = np.atleast_2d(axes) for ax in axes.ravel(): ax.axis("off") for k, (_, c) in enumerate(cands.iterrows()): ax = axes.ravel()[k] im, layer = get_cutout(c.ra, c["dec"], c.label, c[pixscale_col]) if im is not None: ax.imshow(im) h = im.shape[0] # effective radius circle th = np.linspace(0, 2 * np.pi, 60) rpix = c.theta_e / c[pixscale_col] ax.plot(h / 2 + rpix * np.cos(th), h / 2 + rpix * np.sin(th), color="w", lw=0.5, alpha=0.6) cls = c["class"] if "class" in c and isinstance(c["class"], str) else "" mark = "" if "explained" not in c or c.explained == 1 else "UNEXPLAINED: " ax.set_title("%s dI=%.2f sig=%.0f z=%.3f\n%s%s" % (c.label, c.dI, 10 ** c.logs, c.zdist, mark, cls), fontsize=7, fontweight="bold" if mark else "normal") fig.suptitle(title, fontsize=10) fig.tight_layout() fig.savefig(fn, dpi=75) plt.close(fig) def colour_residual(df, opt_mag, w1, extra): """Robust residual of (optical - W1) colour after removing its dependence on size, redshift, sigma.""" col = opt_mag - w1 X = np.column_stack([np.ones(len(df))] + extra) m = np.isfinite(col) & np.all(np.isfinite(X), axis=1) for _ in range(5): c = np.linalg.lstsq(X[m], col[m], rcond=None)[0] res = col - X @ c sd = F.mad_sigma(res[m]) m = np.isfinite(res) & (np.abs(res) < 3 * sd) return res, sd def waste_heat(c, name): """gamma_T = (power of a T-kelvin blackbody normalised to the measured W3 or W4 flux) / (observed bolometric starlight). A galaxy with a fraction f of its starlight absorbed and re-radiated at T would show gamma_T = f/(1-f) (2.3 for f=0.7). Non-detections give upper limits (AllWISE 95% limits). Bolometric corrections (assumed): L_bol = 2.0 nuL_nu(r) or 1.5 nuL_nu(J) for old populations.""" c = c.copy() if name == "sdss": opt = 3631.0 * 10 ** (-0.4 * (c.devmag_r - c.ext_r)) * 4.86e14 * 2.0 else: opt = 1594.0 * 10 ** (-0.4 * c.msel) * 2.43e14 * 1.5 w3 = 31.674 * 10 ** (-0.4 * c.W3mag) * 2.59e13 w4 = 8.363 * 10 ** (-0.4 * c.W4mag) * 1.357e13 c["gamma300K_W3"] = w3 / 0.74 / opt c["gamma150K_W4"] = w4 / 0.74 / opt c["gamma300K_is_limit"] = 1 - c.w3_det c["gamma150K_is_limit"] = 1 - c.w4_det return c def main(): out = {} for name in ("sdss", "6df"): d = pd.read_csv(F.resid_path(name)) w = xmatch(d, "vizier:II/328/allwise", os.path.join(F.RAW, "xmatch_allwise_%s.csv" % name)) w = w.sort_values("angDist").drop_duplicates("id") keep = ["id", "angDist", "W1mag", "W2mag", "W3mag", "W4mag", "e_W3mag", "e_W4mag", "ccf", "ex", "qph"] d = d.merge(w[keep], on="id", how="left") d["w1w2"] = d.W1mag - d.W2mag d["w2w3"] = d.W2mag - d.W3mag d["w3_det"] = d.qph.fillna("UUUU").str[2].isin(list("ABC")).astype(int) d["w4_det"] = d.qph.fillna("UUUU").str[3].isin(list("ABC")).astype(int) if name == "sdss": opt = (d.devmag_r - d.ext_r).values res, sd = colour_residual(d, opt, d.W1mag.values, [np.log10(d.theta_e.values), d.zdist.values, d.logs.values, d.gr.values]) else: opt = d.msel.values res, sd = colour_residual(d, opt, d.W1mag.values, [np.log10(d.theta_e.values), d.zdist.values, d.logs.values]) d["dcol_w1"] = res # >0: optical/NIR flux low relative to W1 (dust, bad optical photometry, W1 excess) d["dcol_w1_nsig"] = res / sd d["label"] = ["%s%d" % ("S" if name == "sdss" else "D", k) for k in range(len(d))] cands = d[d.dI <= CAND_THR].sort_values("dI").reset_index(drop=True) cands["label"] = ["%s%02d" % ("S" if name == "sdss" else "D", k + 1) for k in range(len(cands))] cands["pixscale"] = np.where(cands.theta_e > 12, 1.0, 0.5) if name == "6df" else \ np.where(cands.theta_e > 8, 0.6, 0.3) man = os.path.join(F.HERE, "vetting_manual.csv") if os.path.exists(man): mm = pd.read_csv(man, dtype={"id": str}) mm = mm[mm["sample"] == name][["id", "explained", "class"]] cands = cands.merge(mm, left_on=cands.id.astype(str), right_on="id", how="left", suffixes=("", "_m")) cands = cands.drop(columns=[c for c in ("key_0", "id_m") if c in cands.columns]) # mid-IR waste-heat diagnostic: ratio of 300 K (W3) / 150 K (W4) blackbody power to the observed starlight cands = waste_heat(cands, name) contact_sheet(cands, os.path.join(F.FIGS, "candidates_%s.png" % name), "%s: galaxies with FP surface-brightness deficit dI <= %.2f dex (sorted by dI)" % ("SDSS-MGS own sample" if name == "sdss" else "6dFGSv-J", CAND_THR)) out[name] = (d, cands) frac_w = float(d.W1mag.notna().mean()) print("%s: AllWISE match fraction %.3f; candidates (dI<=%.2f): %d; with W3 detection %d; W2-W3>2.0: %d" % (name, frac_w, CAND_THR, len(cands), cands.w3_det.sum(), ((cands.w2w3 > 2.0) & (cands.w3_det == 1)).sum())) d[["id", "W1mag", "W2mag", "W3mag", "W4mag", "w3_det", "w4_det", "w1w2", "w2w3", "dcol_w1", "dcol_w1_nsig"]] \ .to_csv(os.path.join(F.DATA, "wise_%s.csv" % name), index=False, float_format="%.5g") cols = ["label", "id", "ra", "dec", "zdist", "dI", "dFJ", "logs", "e_logs", "theta_e", "logr", "logi", "msel", "snr", "W1mag", "w1w2", "w2w3", "w3_det", "w4_det", "dcol_w1_nsig"] if name == "sdss": cols += ["logs_port", "ba", "conc", "clean", "ew_ha_port", "d4000", "bpt", "in_howlett", "quality", "plate", "mjd", "fiber", "dmag_dev_petro", "mu50"] else: cols += ["Mt", "Rtd", "sersic_n", "ok_selprob", "ok_chi2", "in_6dFGSv_J", "group"] cols += ["gamma300K_W3", "gamma300K_is_limit", "gamma150K_W4", "gamma150K_is_limit"] cands[cols].to_csv(os.path.join(F.DATA, "candidates_%s_auto.csv" % name), index=False, float_format="%.8g") return out if __name__ == "__main__": sys.exit(0 if main() else 0)