"""Build one-row-per-star tables for each spectroscopic survey (and the Gaia parent sample). Output (pickles in WORK, not in the shared folder; regenerated from RAW): apogee_stars.pkl, lamost_stars.pkl, gspspec_stars.pkl, parent_stars.pkl Every table uses the same column names for Gaia astrometry / photometry, 2MASS and the spectroscopic parameters (teff, logg, feh). No quality cut is applied here except de-duplication; cuts are applied (and counted) in analysis.py. """ import glob import os import sys import numpy as np import pandas as pd from common import RAW, WORK, log, sfd_ebv import download_data as dd PLX_ZP = 0.017 # mas; Gaia EDR3 global parallax zero point (Lindegren et al. 2021), added to parallaxes H_DUST = 100.0 # pc; exponential dust scale height used to scale the SFD column to the star # extinction coefficients per unit E(B-V)_SFD (Schlafly & Finkbeiner 2011 recalibration folded in: x0.86) R_G, R_K, R_BPRP, R_J, R_H = 2.13, 0.306, 1.13, 0.709, 0.449 XM_RENAME = {"Source": "source_id", "RAdeg": "ra", "DEdeg": "dec", "Plx": "plx", "e_Plx": "e_plx", "pmRA": "pmra", "pmDE": "pmde", "RUWE": "ruwe", "epsi": "epsi", "sepsi": "sepsi", "Solved": "solved", "IPDgofha": "ipd_gofha", "IPDfmp": "ipd_fmp", "IPDfow": "ipd_fow", "Dup": "dup", "Gmag": "G", "e_Gmag": "e_G", "BPmag": "BP", "RPmag": "RP", "E(BP/RP)": "bprp_excess", "NBPblend": "nbp_blend", "NRPblend": "nrp_blend", "RV": "rv_gaia", "e_RV": "e_rv_gaia", "o_RV": "n_rv", "RVamp": "rvamp", "vbroad": "vbroad", "GRVSmag": "grvs", "RVSS/N": "rvs_sn", "VarFlag": "varflag", "Teff": "teff_gphot", "logg": "logg_gphot", "AG": "ag_gphot", "RAJ2000": "ra2000", "DEJ2000": "dec2000"} AR_RENAME = {"parallax": "plx", "parallax_error": "e_plx", "pmdec": "pmde", "astrometric_excess_noise": "epsi", "astrometric_excess_noise_sig": "sepsi", "astrometric_params_solved": "solved", "ipd_gof_harmonic_amplitude": "ipd_gofha", "ipd_frac_multi_peak": "ipd_fmp", "ipd_frac_odd_win": "ipd_fow", "duplicated_source": "dup", "phot_g_mean_mag": "G", "phot_bp_mean_mag": "BP", "phot_rp_mean_mag": "RP", "phot_bp_rp_excess_factor": "bprp_excess", "phot_bp_n_blended_transits": "nbp_blend", "phot_rp_n_blended_transits": "nrp_blend", "radial_velocity": "rv_gaia", "radial_velocity_error": "e_rv_gaia", "rv_nb_transits": "n_rv", "rv_amplitude_robust": "rvamp", "grvs_mag": "grvs", "rvs_spec_sig_to_noise": "rvs_sn", "phot_variable_flag": "varflag", "teff_gspphot": "teff_gphot", "logg_gspphot": "logg_gphot", "ag_gspphot": "ag_gphot", "l": "glon", "b": "glat"} def galactic(ra, dec): from astropy.coordinates import SkyCoord import astropy.units as u c = SkyCoord(ra=np.asarray(ra) * u.deg, dec=np.asarray(dec) * u.deg, frame="icrs").galactic return c.l.deg, c.b.deg def add_derived(d): """Distance, extinction estimate, absolute magnitudes and colours (in place).""" if "glon" not in d: d["glon"], d["glat"] = galactic(d.ra.values, d.dec.values) d["plxc"] = d.plx + PLX_ZP d["dist"] = 1000.0 / d.plxc d["sfd"] = sfd_ebv(d.glon.values, d.glat.values) zabs = d.dist * np.abs(np.sin(np.radians(d.glat))) d["dustfrac"] = 1.0 - np.exp(-zabs / H_DUST) d["ebv"] = d.sfd * d.dustfrac # E(B-V)_SFD scaled to the star dm = 5.0 * np.log10(d.dist) - 5.0 d["MG"] = d.G - dm - R_G * d.ebv d["MK"] = d.K - dm - R_K * d.ebv d["bprp0"] = d.BP - d.RP - R_BPRP * d.ebv d["gk0"] = (d.G - d.K) - (R_G - R_K) * d.ebv d["poe"] = d.plx / d.e_plx d["hpx5"] = (d.source_id.values // 34359738368) // (4 ** 7) # HEALPix nested level 5 from source_id return d def attach_2mass(d, name): t = pd.read_csv(os.path.join(WORK, "%s_2mass.csv.gz" % name), low_memory=False) idc = "Source" if "Source" in t.columns else "source_id" t = t.rename(columns={idc: "source_id", "angDist": "d2m", "Jmag": "J", "Hmag": "H", "Kmag": "K", "e_Jmag": "e_J", "e_Hmag": "e_H", "e_Kmag": "e_K", "Qfl": "qfl", "2MASS": "tmass"}) t = t[["source_id", "d2m", "tmass", "J", "H", "K", "e_J", "e_H", "e_K", "qfl", "n2m"]] return d.merge(t, on="source_id", how="left") # ------------------------------------------------------------------ APOGEE STARFLAG_BAD = sum(1 << b for b in (0, 3, 4, 9, 16, 17, 18, 19, 21, 22)) ASPCAP_BAD = sum(1 << b for b in (16, 17, 19, 23, 31)) def build_apogee(): out = os.path.join(WORK, "apogee_stars.pkl") a = pd.concat([pd.read_csv(p, dtype={"GaiaEDR3": str}, low_memory=False) for p in dd.get_apogee()], ignore_index=True) n0 = len(a) a = a[a.GaiaEDR3.notna() & (a.GaiaEDR3.str.strip() != "")].copy() a["source_id"] = a.GaiaEDR3.str.strip().astype(np.int64) sf = a.SFlag.fillna(0).astype(np.int64).values af = a.AFlag.fillna(0).astype(np.int64).values a["sflag"], a["aflag"] = sf, af a["flag_ok"] = (np.bitwise_and(sf, STARFLAG_BAD) == 0) & (np.bitwise_and(af, ASPCAP_BAD) == 0) t1 = a.Ap1T1.fillna(0).astype(np.int64).values t2 = a.Ap2T1.fillna(0).astype(np.int64).values a["main_survey"] = (np.bitwise_and(np.right_shift(t1, 11), 7) > 0) | (np.bitwise_and(np.right_shift(t2, 11), 7) > 0) a["nentry"] = a.groupby("source_id")["source_id"].transform("size") g = a.groupby("source_id") a["teff_spread"] = g["Teff"].transform("max") - g["Teff"].transform("min") a["logg_spread"] = g["logg"].transform("max") - g["logg"].transform("min") a = a.sort_values(["flag_ok", "SNR"], ascending=[False, False]).drop_duplicates("source_id") a = a.rename(columns={"APOGEE": "apogee_id", "Teff": "teff", "e_Teff": "e_teff", "e_logg": "e_logg", "mh": "feh_mh", "feh": "feh_fe", "SNR": "snr", "Jmag": "J", "Hmag": "H", "Ksmag": "K", "e_Jmag": "e_J", "e_Hmag": "e_H", "e_Ksmag": "e_K", "s_HRV": "vscatter", "Nvis": "nvis", "Ncomp": "ncomp", "HRV": "rv_spec"}) a["feh"] = a["feh_mh"] # ASPCAP [M/H] is the metallicity used throughout a["apogee_id"] = a.apogee_id.str.strip() keep = ["source_id", "apogee_id", "Tel", "Field", "teff", "e_teff", "logg", "e_logg", "feh", "e_mh", "feh_fe", "am", "snr", "nvis", "vscatter", "ncomp", "rv_spec", "Vsini", "TeffSp", "loggSp", "sflag", "aflag", "flag_ok", "main_survey", "nentry", "teff_spread", "logg_spread", "Hmin", "Hmax", "jkmin", "jkmax", "J", "H", "K", "e_J", "e_H", "e_K", "AkT", "Survey", "Prog", "EFlag", "Chi2"] a = a[keep] g = pd.read_csv(os.path.join(WORK, "apogee_gaia.csv.gz"), low_memory=False).rename(columns=XM_RENAME) g = g[[c for c in XM_RENAME.values() if c in g.columns]] d = a.merge(g, on="source_id", how="inner") # 2MASS photometric-quality flag, match distance and match multiplicity from the same CDS XMatch as for # the other surveys. The J, H, Ks magnitudes stay those of the APOGEE catalogue (they are 2MASS PSC # values and are what APOGEE targeted on); a star without a 2MASS match within 1.5 arcsec of its Gaia # epoch-2000 position gets no flag and fails the Ks quality cut. t = pd.read_csv(os.path.join(WORK, "apogee_2mass.csv.gz"), low_memory=False) t = t.rename(columns={"Source": "source_id", "angDist": "d2m", "Qfl": "qfl", "2MASS": "tmass", "Kmag": "K_xm"}) d = d.merge(t[["source_id", "d2m", "tmass", "qfl", "n2m", "K_xm"]], on="source_id", how="left") d = add_derived(d) d.to_pickle(out) log("APOGEE: %d catalogue rows -> %d unique Gaia ids -> %d with Gaia DR3 columns" % (n0, len(a), len(d))) return d # ------------------------------------------------------------------ LAMOST def build_lamost(): out = os.path.join(WORK, "lamost_stars.pkl") parts = [pd.read_csv(p, dtype={"GaiaDR3": str}, low_memory=False) for p in dd.get_lamost()] a = pd.concat(parts, ignore_index=True) n0 = len(a) a = a[a.GaiaDR3.notna() & (a.GaiaDR3.str.strip() != "")].copy() a["source_id"] = a.GaiaDR3.str.strip().astype(np.int64) a = a.drop(columns=["GaiaDR3"]) g = a.groupby("source_id") a["nspec"] = g["ObsID"].transform("size") a["teff_med"] = g["teff"].transform("median") a["logg_med"] = g["logg"].transform("median") a["feh_med"] = g["feh"].transform("median") a["teff_spread"] = g["teff"].transform("max") - g["teff"].transform("min") a["logg_spread"] = g["logg"].transform("max") - g["logg"].transform("min") a["rv_spread"] = g["rv_l"].transform("max") - g["rv_l"].transform("min") a = a.sort_values("snrg", ascending=False).drop_duplicates("source_id") a = a.rename(columns={"snrg": "snr", "ObsID": "obsid", "rv_l": "rv_spec"}) gcat = pd.read_csv(os.path.join(WORK, "lamost_gaia.csv.gz"), low_memory=False) gcat["fib_sep"] = gcat["angDist"] # LAMOST fibre position vs Gaia (epoch 2016) position gcat = gcat.rename(columns=XM_RENAME) gcat = gcat[[c for c in list(XM_RENAME.values()) + ["fib_sep"] if c in gcat.columns]] d = a.merge(gcat, on="source_id", how="inner") d = attach_2mass(d, "lamost") d = add_derived(d) d.to_pickle(out) log("LAMOST: %d spectra -> %d unique Gaia ids -> %d with Gaia DR3 columns" % (n0, len(a), len(d))) return d # ------------------------------------------------------------------ Gaia GSP-Spec and parent def _read_gaia(parts): d = pd.concat([pd.read_csv(p, low_memory=False, dtype={"flags_gspspec": str}) for p in parts], ignore_index=True) d = d.rename(columns=AR_RENAME).drop_duplicates("source_id") for c in ("dup", "non_single_star"): if c in d and d[c].dtype == object: d[c] = d[c].astype(str).str.lower().isin(["true", "1", "t"]).astype(int) return d def build_gspspec(): out = os.path.join(WORK, "gspspec_stars.pkl") d = _read_gaia(dd.get_gspspec()) n0 = len(d) d = d.rename(columns={"teff_gspspec": "teff", "logg_gspspec": "logg", "mh_gspspec": "feh"}) d["e_teff"] = 0.5 * (d.teff_gspspec_upper - d.teff_gspspec_lower) d["e_logg"] = 0.5 * (d.logg_gspspec_upper - d.logg_gspspec_lower) d["e_feh"] = 0.5 * (d.mh_gspspec_upper - d.mh_gspspec_lower) d["snr"] = d.rvs_sn d = attach_2mass(d, "gspspec") d = add_derived(d) d.to_pickle(out) log("GSP-Spec: %d rows -> %d stars" % (n0, len(d))) return d def build_parent(): out = os.path.join(WORK, "parent_stars.pkl") d = _read_gaia(dd.get_parent()) d = attach_2mass(d, "parent") d = add_derived(d) d.to_pickle(out) log("Parent Gaia random sample: %d stars" % len(d)) return d BUILDERS = dict(apogee=build_apogee, lamost=build_lamost, gspspec=build_gspspec, parent=build_parent) if __name__ == "__main__": for s in (sys.argv[1:] or list(BUILDERS)): BUILDERS[s]()