"""Download every raw input of the dim_stars study into RAW (idempotent). Sources * APOGEE DR17 allStar VizieR III/286/catalog (TAP) * LAMOST DR11 LRS AFGK stellar VizieR V/162/dr11sl (TAP) * Gaia DR3 GSP-Spec + astrometry Gaia archive, gaiadr3.astrophysical_parameters x gaia_source (synchronous TAP, chunked by HEALPix ranges of source_id) * Gaia DR3 random parent sample Gaia archive, 2% of gaia_source by random_index (synchronous TAP, chunked) * Gaia DR3 columns for APOGEE/LAMOST stars CDS XMatch vs I/355/gaiadr3 * 2MASS PSC CDS XMatch vs II/246/out (all four samples, APOGEE included) * SFD E(B-V) maps github.com/kbarbary/sfddata Run: python3 download_data.py [sfd|apogee|lamost|gspspec|parent|xmatch_apogee|xmatch_lamost|xmatch_gspspec|xmatch_parent ...] """ import os import sys import numpy as np import pandas as pd from common import (RAW, WORK, VIZIER_TAP, GAIA_TAP, log, tap_sync_to_file, gaia_async_to_file, xmatch_upload, sfd_ebv, to_epoch2000) TEFF_LO, TEFF_HI, LOGG_LO = 3900, 7100, 3.4 GAIA_XCOLS = ("Source,RAdeg,DEdeg,RAJ2000,DEJ2000,Plx,e_Plx,pmRA,pmDE,RUWE,epsi,sepsi,Solved,IPDgofha,IPDfmp," "IPDfow,Dup,Gmag,e_Gmag,BPmag,RPmag,E(BP/RP),NBPblend,NRPblend,RV,e_RV,o_RV,RVamp,vbroad," "GRVSmag,RVSS/N,VarFlag,Teff,logg,AG") # ---------------------------------------------------------------- APOGEE DR17 def get_apogee(): cols = ('APOGEE, Tel, Field, RAJ2000, DEJ2000, GLON, GLAT, Jmag, e_Jmag, Hmag, e_Hmag, Ksmag, e_Ksmag, ' 'AkT, "E(B-V)" AS sfd_ebv, Ap1T1, Ap1T2, Ap2T1, Ap2T2, Ap2T3, Survey, Prog, Nvis, SNR, SFlag, SaFlag, ' 'HRV, s_HRV, e_HRV, Ncomp, Hmin, Hmax, "J-Kmin" AS jkmin, "J-Kmax" AS jkmax, GaiaEDR3, plx, e_plx, ' 'Gmag, BPmag, RPmag, Chi2, AFlag, EFlag, Teff, e_Teff, logg, e_logg, "[M/H]" AS mh, "e_[M/H]" AS e_mh, ' '"[a/M]" AS am, Vmicro, Vsini, TeffSp, loggSp, "[Fe/H]" AS feh, "f_[Fe/H]" AS f_feh') parts = [] edges = [0, 60, 120, 180, 240, 300, 360.0001] for a, b in zip(edges[:-1], edges[1:]): p = os.path.join(RAW, "apogee_dr17_ra%03d.csv" % a) q = ('SELECT %s FROM "III/286/catalog" WHERE logg >= %g AND Teff BETWEEN %d AND %d ' 'AND RAJ2000 >= %g AND RAJ2000 < %g' % (cols, LOGG_LO, TEFF_LO, TEFF_HI, a, b)) log("APOGEE chunk RA", a) tap_sync_to_file(VIZIER_TAP, q, p, maxrec=2000000) parts.append(p) return parts # ---------------------------------------------------------------- LAMOST DR11 def get_lamost(step=10): cols = ('l.ObsID, l.RAJ2000 AS ral, l.DEJ2000 AS decl, l.snrg, l.snrr, l.snri, l.subClass, l.Teff AS teff, ' 'l.e_Teff AS e_teff, l.logg AS logg, l.e_logg AS e_logg, l."[Fe/H]" AS feh, l."e_[Fe/H]" AS e_feh, ' 'l.RV AS rv_l, l.e_RV AS e_rv_l, l.offsets, l.GaiaDR3') parts = [] for a in range(0, 360, step): p = os.path.join(RAW, "lamost_dr11_ra%03d.csv" % a) q = ('SELECT %s FROM "V/162/dr11sl" AS l WHERE l.RAJ2000 >= %d AND l.RAJ2000 < %d AND l.snrg >= 30 ' 'AND l.logg >= %g AND l.Teff BETWEEN %d AND %d' % (cols, a, a + step, LOGG_LO, TEFF_LO, TEFF_HI)) if not os.path.exists(p): log("LAMOST chunk RA", a) tap_sync_to_file(VIZIER_TAP, q, p, maxrec=3000000) parts.append(p) return parts # ---------------------------------------------------------------- Gaia GSP-Spec GAIA_COLS = ("g.source_id, g.ra, g.dec, g.l, g.b, g.parallax, g.parallax_error, g.pmra, g.pmdec, g.ruwe, " "g.astrometric_excess_noise, g.astrometric_excess_noise_sig, g.astrometric_params_solved, " "g.ipd_gof_harmonic_amplitude, g.ipd_frac_multi_peak, g.ipd_frac_odd_win, g.duplicated_source, " "g.phot_g_mean_mag, g.phot_g_mean_flux_over_error, g.phot_bp_mean_mag, g.phot_rp_mean_mag, " "g.phot_bp_rp_excess_factor, g.phot_bp_n_blended_transits, g.phot_rp_n_blended_transits, " "g.radial_velocity, g.radial_velocity_error, g.rv_nb_transits, g.rv_amplitude_robust, g.vbroad, " "g.grvs_mag, g.rvs_spec_sig_to_noise, g.phot_variable_flag, g.non_single_star, " "g.teff_gspphot, g.logg_gspphot, g.ag_gspphot") def _gaia_sync_chunks(name, queries, nrow_warn=1500000): """Run a list of (label, adql) synchronous Gaia-archive queries and concatenate into one CSV. (Asynchronous anonymous jobs were queueing for many minutes when this study was run; small synchronous queries driven by an indexed column return in seconds.)""" out = os.path.join(RAW, name) if os.path.exists(out) and os.path.getsize(out) > 0: return out tmpdir = os.path.join(RAW, name + ".parts") os.makedirs(tmpdir, exist_ok=True) files = [] for lab, q in queries: p = os.path.join(tmpdir, "%s.csv" % lab) if not os.path.exists(p): tap_sync_to_file(GAIA_TAP + "/sync", q, p, maxrec=3000000, timeout=600) files.append(p) with open(out + ".part", "w") as fo: first = True for p in files: with open(p) as fi: head = fi.readline() if first: fo.write(head) first = False n = 0 for line in fi: fo.write(line) n += 1 if n >= nrow_warn: raise RuntimeError("chunk %s may be truncated (%d rows)" % (p, n)) os.replace(out + ".part", out) for p in files: # the chunks are now redundant (about 1 GB in total) os.remove(p) try: os.rmdir(tmpdir) except OSError: pass return out def get_gspspec(): """All GSP-Spec FGK dwarfs/subgiants with parallax > 1 mas. Split by HEALPix level-2 pixel (source_id range) so that every query is a short index range scan of astrophysical_parameters.""" parts = [] pixels = [int(x) for x in os.environ.get("GSP_PIX", ",".join(map(str, range(12)))).split(",")] for pix in pixels: qs = [] for sub in range(16): lo = (pix * 16 + sub) * 2 ** 55 hi = lo + 2 ** 55 - 1 q = ("SELECT " + GAIA_COLS + ", a.teff_gspspec, a.teff_gspspec_lower, a.teff_gspspec_upper, " "a.logg_gspspec, a.logg_gspspec_lower, a.logg_gspspec_upper, a.mh_gspspec, a.mh_gspspec_lower, " "a.mh_gspspec_upper, a.alphafe_gspspec, a.logchisq_gspspec, a.flags_gspspec " "FROM gaiadr3.astrophysical_parameters AS a JOIN gaiadr3.gaia_source AS g USING (source_id) " "WHERE a.source_id BETWEEN %d AND %d AND a.teff_gspspec >= %d AND a.teff_gspspec <= %d " "AND a.logg_gspspec >= %g AND g.parallax > 1.0 AND g.parallax_over_error > 10" % (lo, hi, TEFF_LO, TEFF_HI, LOGG_LO)) qs.append(("sub%02d" % sub, q)) if not os.path.exists(os.path.join(RAW, "gaia_gspspec_hpx%02d.csv" % pix)): log("GSP-Spec HEALPix base pixel", pix) parts.append(_gaia_sync_chunks("gaia_gspspec_hpx%02d.csv" % pix, qs)) return parts def get_parent(): """2% random sample of Gaia DR3 (random_index), used only for the selection-function model.""" parts = [] n = 1811709771 step = int(0.005 * n) for k in range(4): qs = [] nsub = 6 for j in range(nsub): lo = k * step + j * (step // nsub) hi = k * step + (j + 1) * (step // nsub) if j < nsub - 1 else (k + 1) * step q = ("SELECT " + GAIA_COLS + " FROM gaiadr3.gaia_source AS g WHERE g.random_index >= %d AND " "g.random_index < %d AND g.parallax > 1.0 AND g.parallax_over_error > 10 AND " "g.phot_g_mean_mag < 17.5" % (lo, hi)) qs.append(("sub%02d" % j, q)) if not os.path.exists(os.path.join(RAW, "gaia_parent_%d.csv" % k)): log("parent slice", k) parts.append(_gaia_sync_chunks("gaia_parent_%d.csv" % k, qs)) return parts # ---------------------------------------------------------------- cross-matches def read_ids(path, usecols, idcol): """Read a CSV keeping 19-digit Gaia source ids exact (never via float64).""" d = pd.read_csv(path, usecols=usecols, dtype={idcol: str}) d = d[d[idcol].notna() & (d[idcol].str.strip() != "")].copy() d[idcol] = d[idcol].str.strip().astype(np.int64) return d def _chunks(df, n): for i in range(0, len(df), n): yield df.iloc[i:i + n] def xmatch_gaia_for(name, df, idcol, racol, deccol, radius=5.0, chunk=250000): """Fetch Gaia DR3 columns for stars with a known Gaia source_id by positional match + id equality.""" out = os.path.join(WORK, "%s_gaia.parquet" % name) out = out.replace(".parquet", ".csv.gz") if os.path.exists(out): return out res = [] for i, c in enumerate(_chunks(df[[idcol, racol, deccol]], chunk)): log(" xmatch Gaia %s chunk %d (%d rows)" % (name, i, len(c))) x = xmatch_upload(c, racol, deccol, "vizier:I/355/gaiadr3", radius, cols2=GAIA_XCOLS) x = x[x["Source"] == x[idcol]].drop_duplicates("Source") res.append(x) x = pd.concat(res, ignore_index=True).drop_duplicates("Source") x.to_csv(out + ".tmp.csv.gz", index=False) os.replace(out + ".tmp.csv.gz", out) log(" -> %s rows with Gaia DR3 columns (%d requested)" % (len(x), len(df))) return out def xmatch_2mass_for(name, df, idcol, racol, deccol, radius=1.5, chunk=400000): """2MASS PSC match using positions at epoch 2000 (best = nearest; number of matches kept).""" out = os.path.join(WORK, "%s_2mass.csv.gz" % name) if os.path.exists(out): return out up = df[[idcol, racol, deccol]].rename(columns={racol: "ra_up", deccol: "dec_up"}) res = [] for i, c in enumerate(_chunks(up, chunk)): log(" xmatch 2MASS %s chunk %d (%d rows)" % (name, i, len(c))) x = xmatch_upload(c, "ra_up", "dec_up", "vizier:II/246/out", radius, cols2="2MASS,Jmag,Hmag,Kmag,e_Jmag,e_Hmag,e_Kmag,Qfl,Rfl,X") x["n2m"] = x.groupby(idcol)[idcol].transform("size") x = x.sort_values("angDist").drop_duplicates(idcol) res.append(x.drop(columns=["ra_up", "dec_up"])) x = pd.concat(res, ignore_index=True).drop_duplicates(idcol) x.to_csv(out + ".tmp.csv.gz", index=False) os.replace(out + ".tmp.csv.gz", out) log(" -> %d 2MASS matches (%d requested)" % (len(x), len(df))) return out def xmatch_lamost(): lam = pd.concat([read_ids(p, ["ral", "decl", "GaiaDR3"], "GaiaDR3") for p in get_lamost()], ignore_index=True) lamu = lam.drop_duplicates("GaiaDR3") log("LAMOST spectra %d, unique Gaia ids %d" % (len(lam), len(lamu))) g = pd.read_csv(xmatch_gaia_for("lamost", lamu, "GaiaDR3", "ral", "decl"), low_memory=False) xmatch_2mass_for("lamost", g[["Source", "RAJ2000", "DEJ2000"]], "Source", "RAJ2000", "DEJ2000") def xmatch_apogee(): apo = pd.concat([read_ids(p, ["RAJ2000", "DEJ2000", "GaiaEDR3"], "GaiaEDR3") for p in get_apogee()], ignore_index=True) apo = apo.rename(columns={"RAJ2000": "ra_ap", "DEJ2000": "dec_ap"}).drop_duplicates("GaiaEDR3") g = pd.read_csv(xmatch_gaia_for("apogee", apo, "GaiaEDR3", "ra_ap", "dec_ap"), low_memory=False) # 2MASS quality flags and match multiplicity, exactly as for the other surveys (the APOGEE catalogue # carries 2MASS magnitudes but not the photometric-quality flag) xmatch_2mass_for("apogee", g[["Source", "RAJ2000", "DEJ2000"]], "Source", "RAJ2000", "DEJ2000") def xmatch_gaia_samples(which=("gspspec", "parent")): # GSP-Spec and parent: 2MASS via epoch-2000 positions getters = dict(gspspec=get_gspspec, parent=get_parent) for name in which: parts = getters[name]() d = pd.concat([pd.read_csv(p, usecols=["source_id", "ra", "dec", "pmra", "pmdec"]) for p in parts], ignore_index=True).drop_duplicates("source_id") d["ra2000"], d["dec2000"] = to_epoch2000(d.ra.values, d.dec.values, d.pmra.values, d.pmdec.values) xmatch_2mass_for(name, d[["source_id", "ra2000", "dec2000"]], "source_id", "ra2000", "dec2000") def main(steps): if not steps: steps = ["sfd", "apogee", "lamost", "gspspec", "parent", "xmatch_apogee", "xmatch_lamost", "xmatch_gspspec", "xmatch_parent"] for s in steps: log("=== download step:", s) if s == "sfd": log("SFD test E(B-V) at (l,b)=(0,90),(120,-30):", sfd_ebv([0.0, 120.0], [90.0, -30.0])) elif s == "apogee": get_apogee() elif s == "lamost": get_lamost() elif s == "gspspec": get_gspspec() elif s == "parent": get_parent() elif s == "xmatch_apogee": xmatch_apogee() elif s == "xmatch_lamost": xmatch_lamost() elif s == "xmatch_gspspec": xmatch_gaia_samples(("gspspec",)) elif s == "xmatch_parent": xmatch_gaia_samples(("parent",)) else: raise SystemExit("unknown step " + s) log("downloads complete") if __name__ == "__main__": main(sys.argv[1:])