"""Step 6a: external-catalogue queries for every excess star (cached in RAW/vet/). For each detected excess star we collect * SIMBAD identification through its Gaia DR3 identifier: main type, all object types, spectral type, parent associations (h_link); * everything within 12-20 arcsec of the position WISE saw it at (epoch 2010.5) in SIMBAD, Gaia DR3, Pan-STARRS1, SDSS DR16, DES DR2, CatWISE2020, 2MASS XSC and the radio surveys NVSS, FIRST, VLASS-QL, RACS-low and SUMSS (CDS XMatch service); * membership of published infrared-excess / debris-disc / young-association catalogues (5 arcsec around the epoch-2000 position). """ import io import os import sys import numpy as np import pandas as pd from common import RAW, DATA, SIMBAD_TAP, tap_sync, cds_xmatch import importlib inj = importlib.import_module("04_injection") VET = os.path.join(RAW, "vet") VET_CMP = os.path.join(RAW, "vet_cmp") N_CMP = 10 # comparison (non-excess) stars per excess star XM_WISE_EPOCH = { # name: (cat2, radius arcsec) "simbad": ("simbad", 12), "gaia": ("vizier:I/355/gaiadr3", 12), "ps1": ("vizier:II/349/ps1", 10), "sdss": ("vizier:V/154/sdss16", 10), "des": ("vizier:II/371/des_dr2", 10), "catwise": ("vizier:II/365/catwise", 12), "xsc": ("vizier:VII/233/xsc", 20), "nvss": ("vizier:VIII/65/nvss", 15), "first": ("vizier:VIII/92/first14", 10), "vlass": ("vizier:J/ApJS/255/30/comp", 10), "racs": ("vizier:J/other/PASA/38.58/gausscut", 12), "sumss": ("vizier:VIII/81B/sumss212", 15), } XM_J2000 = { "cotten_prime": ("vizier:J/ApJS/225/15/table3", 5), "cotten_reserved": ("vizier:J/ApJS/225/15/table4", 5), "patel14": ("vizier:J/ApJS/212/10/catalog", 5), "chen14": ("vizier:J/ApJS/211/25/catalog", 5), "mcdonald17": ("vizier:J/MNRAS/471/770/table3", 5), "theissen14": ("vizier:J/ApJ/794/146/catalog", 5), "theissen17": ("vizier:J/AJ/153/165/table4", 5), "gagne18_members": ("vizier:J/ApJ/856/23/members", 5), "gagne18b_candidates": ("vizier:J/ApJ/862/138/table2", 5), "gagne18b_bonafide": ("vizier:J/ApJ/862/138/table5", 5), } def _write_ids(path, ids): with open(path, "w") as f: f.write("\n".join(str(int(i)) for i in ids) + "\n") def _read_ids(path): with open(path) as f: return {int(x) for x in f.read().split()} def positions(c): """Positions at epochs 2000.0, 2010.5 (WISE) and 2016.0 (Gaia) from Gaia astrometry.""" out = pd.DataFrame({"sid": c.sid.values}) cosd = np.cos(np.radians(c.DE_ICRS.values)) for name, ep in (("2000", 2000.0), ("wise", 2010.5), ("2016", 2016.0)): dt = ep - 2016.0 out["ra_" + name] = (c.RA_ICRS.values + c.pmRA.fillna(0).values * dt / 3.6e6 / cosd) % 360.0 out["dec_" + name] = c.DE_ICRS.values + c.pmDE.fillna(0).values * dt / 3.6e6 return out def run_xmatch(pos, tag, table, name, force=False, vet_dir=None): VET = vet_dir or globals()["VET"] os.makedirs(VET, exist_ok=True) out = os.path.join(VET, f"xm_{name}.csv") if os.path.exists(out) and not force: donef = os.path.join(VET, f"xm_{name}.sids") if os.path.exists(donef) and _read_ids(donef) >= set(int(x) for x in pos.sid): return cat2, rad = table up = pos[["sid", "ra_" + tag, "dec_" + tag]].rename(columns={"ra_" + tag: "ra", "dec_" + tag: "dec"}) txt = cds_xmatch(up.to_csv(index=False, float_format="%.7f"), cat2, rad) if not txt.startswith("angDist"): raise RuntimeError(f"XMatch {name} failed: {txt[:300]}") with open(out, "w") as f: f.write(txt) _write_ids(os.path.join(VET, f"xm_{name}.sids"), pos.sid.values) print(f"xmatch {name}: {txt.count(chr(10)) - 1} rows", flush=True) def simbad_by_id(sids, force=False): out = os.path.join(VET, "simbad_basic.csv") donef = os.path.join(VET, "simbad_basic.sids") if os.path.exists(out) and os.path.exists(donef) and not force: if _read_ids(donef) >= set(int(x) for x in sids): return basic, otypes, parents = [], [], [] sids = list(sids) for k in range(0, len(sids), 100): ids = ",".join(f"'Gaia DR3 {s}'" for s in sids[k:k + 100]) basic.append(pd.read_csv(io.StringIO(tap_sync( SIMBAD_TAP, "SELECT i.id AS gaia_id, b.oid, b.main_id, b.otype, b.sp_type, b.nbref, b.ra, b.dec " f"FROM ident AS i JOIN basic AS b ON b.oid = i.oidref WHERE i.id IN ({ids})")))) otypes.append(pd.read_csv(io.StringIO(tap_sync( SIMBAD_TAP, f"SELECT i.id AS gaia_id, o.otype FROM ident AS i JOIN otypes AS o ON o.oidref = i.oidref " f"WHERE i.id IN ({ids})")))) parents.append(pd.read_csv(io.StringIO(tap_sync( SIMBAD_TAP, "SELECT i.id AS gaia_id, p.main_id AS parent, p.otype AS parent_otype, h.membership " "FROM ident AS i JOIN h_link AS h ON h.child = i.oidref JOIN basic AS p ON p.oid = h.parent " f"WHERE i.id IN ({ids})")))) for name, lst in (("simbad_basic", basic), ("simbad_otypes", otypes), ("simbad_parents", parents)): df = pd.concat(lst, ignore_index=True) df["sid"] = df.gaia_id.str.replace("Gaia DR3 ", "").astype("int64") df.drop_duplicates().to_csv(os.path.join(VET, name + ".csv"), index=False) _write_ids(donef, sids) print("simbad by id done", flush=True) SGRO_BIBCODE = "2021MNRAS.508.3084S" # Sgro & Song 2021: W3/W4 excess among Gaia DR2 K/M dwarfs within 100 pc def simbad_by_bibcode(bibcode=SGRO_BIBCODE, name="simbad_sgro2021", force=False): """Objects that SIMBAD attaches to a paper, with their Gaia DR3 identifier. Sgro & Song (2021) have no VizieR table, so this is the only machine-readable list of the stars they discuss (their 37 new candidates, the 12 known excess stars they re-analyse, and comparison objects).""" out = os.path.join(VET, name + ".csv") if os.path.exists(out) and os.path.getsize(out) > 50 and not force: return txt = tap_sync(SIMBAD_TAP, "SELECT b.main_id, b.otype, b.ra, b.dec, i.id AS gaia_id FROM ref AS r " "JOIN has_ref AS h ON h.oidbibref = r.oidbib JOIN basic AS b ON b.oid = h.oidref " "JOIN ident AS i ON i.oidref = b.oid " f"WHERE r.bibcode = '{bibcode}' AND i.id LIKE 'Gaia DR3 %'") df = pd.read_csv(io.StringIO(txt)) df["sid"] = df.gaia_id.str.replace("Gaia DR3 ", "").astype("int64") df.to_csv(out, index=False) print(f"SIMBAD objects of {bibcode}: {len(df)}", flush=True) def comparison_set(d, n_per=N_CMP): """For every excess star: the n_per nearest (on the sky) searchable NON-excess stars of the same class and similar brightness (|dG| < 1.5). Used to measure how often the background criteria of step 7 fire on stars with no excess (base rate of chance association).""" from scipy.spatial import cKDTree def unit(ra, dec): ra, dec = np.radians(ra), np.radians(dec) return np.c_[np.cos(dec) * np.cos(ra), np.cos(dec) * np.sin(ra), np.sin(dec)] det = d[d.det] rows, used = [], set() for cls in ("A", "FGK", "M"): pool = d[(d.cls == cls) & (d.q3 | d.q4) & ~d.det & (d.chi3.abs().fillna(0) < 3) & (d.chi4.abs().fillna(0) < 3)] tree = cKDTree(unit(pool.RA_ICRS.values, pool.DE_ICRS.values)) psid, pg = pool.sid.values, pool.Gmag.values for _, s in det[det.cls == cls].sort_values("sid").iterrows(): _, idx = tree.query(unit(np.array([s.RA_ICRS]), np.array([s.DE_ICRS]))[0], k=min(400, len(pool))) got = 0 for j in np.atleast_1d(idx): if psid[j] in used or abs(pg[j] - s.Gmag) > 1.5: continue used.add(psid[j]) rows.append(dict(sid=psid[j], match_sid=s.sid)) got += 1 if got == n_per: break m = pd.DataFrame(rows) return d[d.sid.isin(m.sid)].merge(m, on="sid") def main(force=False): d = inj.load_stars() c = d[d.det].copy() pos = positions(c) simbad_by_id(pos.sid.values, force) simbad_by_bibcode(force=force) for name, table in XM_WISE_EPOCH.items(): run_xmatch(pos, "wise", table, name, force) for name, table in XM_J2000.items(): run_xmatch(pos, "2000", table, name, force) pos.to_csv(os.path.join(VET, "positions.csv"), index=False, float_format="%.7f") print("vetting queries done for", len(pos), "excess stars") # comparison (non-excess) stars: same positional cross-matches cmp_ = comparison_set(d) cpos = positions(cmp_) cpos["match_sid"] = cmp_.match_sid.values for name, table in XM_WISE_EPOCH.items(): run_xmatch(cpos, "wise", table, name, force, vet_dir=VET_CMP) cpos.to_csv(os.path.join(VET_CMP, "positions.csv"), index=False, float_format="%.7f") print("comparison queries done for", len(cpos), "non-excess stars") if __name__ == "__main__": main(force="--force" in sys.argv)