#!/usr/bin/env python3 """Merge the raw downloads into one per-galaxy master table. CACHE/master.csv all columns (large; re-creatable) data/master_slim.csv.gz the columns used by the analysis (kept in the study folder) data/build_summary.json row counts at every step Photometry columns (all flux densities in mJy on the WISE Vega-system scale F0*10^(-0.4 m)): F1..F4, e1..e4 AllWISE profile-fit ("raw" w?flux / w?sigflux; negative values kept) L1..L4, le1..le4 Legacy Surveys DR9 Tractor forced photometry on unWISE coadds (optical-model shapes) G1, G3, G4, ge3, ge4 AllWISE elliptical-aperture (2MASS XSC shape) photometry where measured (w?gmag) """ import json import os import numpy as np import pandas as pd import warnings warnings.simplefilter("ignore", pd.errors.PerformanceWarning) import mirq_common as M # MPA-JHU recommended scalings of the formal line-flux uncertainties (Brinchmann et al. 2013, DR8 notes) ERR_SCALE = dict(h_alpha=2.473, h_beta=1.882, oiii_5007=1.566, nii_6584=2.039) def line_classes(d): """Transparent emission-line classes from MPA-JHU measurements inside the 3-arcsec fibre. nodata no usable MPA-JHU H-alpha measurement: no row, or a failed fit (flux error <= 0; the failed fits have flux = 0 and a negative or zero error). Revision: the first pass called the failed fits 'lineless', which let them into the clean subsample without a line measurement. lineless H-alpha not detected (S/N < 3) or emission EW < 0.5 A retired H-alpha detected, 0.5 <= EW < 3 A (weak LINER-like emission that old stars can power; Cid Fernandes et al. 2011) -- NOT counted as an explanation of a mid-IR excess SF / composite / AGN EW >= 3 A; BPT position (Kauffmann+03 and Kewley+01 lines) if all four lines have S/N >= 3, otherwise [NII]/H-alpha alone (log < -0.4 SF, > -0.2 AGN, between composite) """ sn = {} for k, s in ERR_SCALE.items(): with np.errstate(divide="ignore", invalid="ignore"): sn[k] = d[k + "_flux"] / (s * d[k + "_flux_err"].where(d[k + "_flux_err"] > 0)) ew = -d.h_alpha_eqw # positive = emission cls = pd.Series("lineless", index=d.index) det = (sn["h_alpha"] >= 3) & (ew >= 0.5) cls[det & (ew < 3)] = "retired" strong = det & (ew >= 3) with np.errstate(divide="ignore", invalid="ignore"): n2 = np.log10(d.nii_6584_flux / d.h_alpha_flux) o3 = np.log10(d.oiii_5007_flux / d.h_beta_flux) four = strong & (sn["h_beta"] >= 3) & (sn["oiii_5007"] >= 3) & (sn["nii_6584"] >= 3) kauf = (n2 < 0.05) & (o3 < 0.61 / (n2 - 0.05) + 1.3) kew = (n2 < 0.47) & (o3 < 0.61 / (n2 - 0.47) + 1.19) cls[four & kauf] = "SF" cls[four & ~kauf & kew] = "composite" cls[four & ~kew] = "AGN" two = strong & ~four cls[two & (n2 < -0.4)] = "SF" cls[two & (n2 >= -0.4) & (n2 <= -0.2)] = "composite" cls[two & ((n2 > -0.2) | ~np.isfinite(n2))] = "AGN" cls[d.h_alpha_flux.isna() | ~(d.h_alpha_flux_err > 0)] = "nodata" return cls, sn["h_alpha"], n2, o3 def main(): summ = {} s = pd.read_csv(os.path.join(M.DATA, "sample_etg.csv.gz")) summ["n_sample"] = len(s) # ---------------- AllWISE: nearest source within 6 arcsec ---------------- aw = pd.read_csv(os.path.join(M.RAW, "allwise_match.csv"), dtype={"cc_flags": str, "ph_qual": str}) aw = aw.rename(columns={"sid": "id", "ra": "ra_w", "dec": "dec_w"}) aw = aw.merge(s[["id", "ra", "dec"]], on="id", how="left") aw["sep_w"] = M.angsep_arcsec(aw.ra.values, aw["dec"].values, aw.ra_w.values, aw.dec_w.values) nmatch = aw.groupby("id").size().rename("n_allwise_6as") aw = aw.sort_values("sep_w").drop_duplicates("id").drop(columns=["ra", "dec"]) d = s.merge(aw, on="id", how="left").merge(nmatch, on="id", how="left") d["n_allwise_6as"] = d.n_allwise_6as.fillna(0).astype(int) summ["n_allwise_any_6as"] = int((d.n_allwise_6as > 0).sum()) summ["n_allwise_within_3as"] = int((d.sep_w <= 3).sum()) for b in (1, 2, 3, 4): d["F%d" % b] = M.wise_dn_to_jy(d["w%dflux" % b], b) * 1e3 d["e%d" % b] = M.wise_dn_to_jy(d["w%dsigflux" % b], b) * 1e3 for b in (1, 3, 4): g = M.wise_mag_to_jy(d["w%dgmag" % b], b) * 1e3 d["G%d" % b] = g d["ge%d" % b] = g * d["w%dgerr" % b] * np.log(10) / 2.5 a = M.wise_mag_to_jy(d["w%dmag_8" % b], b) * 1e3 d["A%d" % b] = a d["ae%d" % b] = a * d["w%dsigm_8" % b] * np.log(10) / 2.5 # ---------------- Legacy Surveys DR9 forced photometry ---------------- ls = pd.read_csv(os.path.join(M.RAW, "lsdr9_match.csv")) ls = ls.rename(columns={"specobjid": "specObjID", "distance": "sep_ls", "type": "ls_type"}) for b in (1, 2, 3, 4): ls["L%d" % b] = M.ls_nmgy_to_wise_jy(ls["flux_w%d" % b], b) * 1e3 with np.errstate(divide="ignore"): ls["le%d" % b] = M.ls_nmgy_to_wise_jy(1.0 / np.sqrt(ls["flux_ivar_w%d" % b].where(ls["flux_ivar_w%d" % b] > 0)), b) * 1e3 keep = ["specObjID", "sep_ls", "ls_id", "ls_type", "shape_r", "sersic", "flux_g", "flux_r", "flux_z", "fracflux_w1", "fracflux_w3", "fracflux_w4", "rchisq_w1", "rchisq_w3", "rchisq_w4", "nobs_w3", "nobs_w4", "maskbits", "fracflux_r"] + ["L%d" % b for b in (1, 2, 3, 4)] + ["le%d" % b for b in (1, 2, 3, 4)] d = d.merge(ls[keep], on="specObjID", how="left") summ["n_lsdr9"] = int(d.L1.notna().sum()) # ---------------- SDSS MPA-JHU lines ---------------- ln = pd.read_csv(os.path.join(M.RAW, "sdss_lines.csv")) d = d.merge(ln, on="specObjID", how="left") cls, sn_ha, n2, o3 = line_classes(d) d["line_class"], d["sn_ha"], d["log_n2ha"], d["log_o3hb"] = cls, sn_ha, n2, o3 d["ew_ha"] = -d.h_alpha_eqw summ["line_class_counts"] = cls.value_counts().to_dict() summ["n_halpha_fit_failed"] = int((d.h_alpha_flux.notna() & ~(d.h_alpha_flux_err > 0)).sum()) # ---------------- radio ---------------- fi = pd.read_csv(os.path.join(M.RAW, "xmatch_first.csv")).sort_values("angDist") fcore = fi[fi.angDist <= 3.0].drop_duplicates("id")[["id", "angDist", "Fpeak", "Fint"]] fcore.columns = ["id", "first_sep", "first_fpeak", "first_fint"] fsum = fi.groupby("id").agg(first_n10=("Fint", "size"), first_fint10=("Fint", "sum")).reset_index() nv = pd.read_csv(os.path.join(M.RAW, "xmatch_nvss.csv")).sort_values("angDist") nv = nv[nv.angDist <= 20.0].drop_duplicates("id")[["id", "angDist", "S1.4"]] nv.columns = ["id", "nvss_sep", "nvss_s14"] d = d.merge(fcore, on="id", how="left").merge(fsum, on="id", how="left").merge(nv, on="id", how="left") summ["n_first_3as"] = int(d.first_fint.notna().sum()) summ["n_first_10as"] = int(d.first_n10.notna().sum()) summ["n_nvss_20as"] = int(d.nvss_s14.notna().sum()) # ---------------- GALEX ---------------- gx = pd.read_csv(os.path.join(M.RAW, "xmatch_galex.csv")).sort_values("angDist").drop_duplicates("id") gx = gx.rename(columns={"angDist": "galex_sep", "E(B-V)": "galex_ebv"})[ ["id", "galex_sep", "NUVmag", "e_NUVmag", "FUVmag", "e_FUVmag", "galex_ebv"]] d = d.merge(gx, on="id", how="left") # Galactic extinction: A_NUV = 8.2 E(B-V) (Wyder et al. 2007); r from SDSS d["nuv_r"] = (d.NUVmag - 8.2 * d.galex_ebv) - (d.modelMag_r - d.extinction_r) summ["n_galex_nuv"] = int(d.NUVmag.notna().sum()) d.to_csv(os.path.join(M.CACHE, "master.csv"), index=False) slim = ["id", "ra", "dec", "z", "zhelio", "veldisp", "theta_e", "Mr", "msel", "snr", "ba", "gr", "d4000", "sep_w", "n_allwise_6as", "ext_flg", "cc_flags", "ph_qual", "nb", "w1mpro", "w2mpro", "w3mpro", "w4mpro", "w1snr", "w2snr", "w3snr", "w4snr", "w3cov", "w4cov", "F1", "F2", "F3", "F4", "e1", "e2", "e3", "e4", "G1", "G3", "G4", "ge3", "ge4", "L1", "L2", "L3", "L4", "le1", "le2", "le3", "le4", "sep_ls", "ls_type", "shape_r", "sersic", "fracflux_w3", "fracflux_w4", "maskbits", "flux_r", "line_class", "ew_ha", "sn_ha", "log_n2ha", "log_o3hb", "lick_hd_a", "modelMag_r", "extinction_r", "first_sep", "first_fpeak", "first_fint", "first_n10", "first_fint10", "nvss_sep", "nvss_s14", "NUVmag", "e_NUVmag", "nuv_r"] out = d[slim].copy() for c in ("ra", "dec"): # keep 0.04-arcsec precision on positions, 4 digits elsewhere out[c] = out[c].map(lambda v: "%.5f" % v) out.to_csv(os.path.join(M.DATA, "master_slim.csv.gz"), index=False, float_format="%.4g", compression="gzip") json.dump(summ, open(os.path.join(M.DATA, "build_summary.json"), "w"), indent=1) print(json.dumps(summ, indent=1)) if __name__ == "__main__": main()