#!/usr/bin/env python3 """Download every raw input of the mir_quiescent study into RAW (cached; safe to re-run). 1. sample data/sample_etg.csv.gz: the 114,210 SDSS early-type galaxies of studies/fp_kiii (slim copy). If the file is missing it is rebuilt from the fp_kiii cache or, failing that, from SDSS SkyServer with the fp_kiii code kept verbatim in code/fp_kiii_build/. 2. sdss_lines SDSS DR17 SkyServer: MPA-JHU line fluxes / classes / indices + ugriz model magnitudes 3. allwise IRSA TAP (allwise_p3as_psd) positional upload match, 6 arcsec, all WISE photometry flavours 4. lsdr9 NOIRLab Astro Data Lab TAP: Legacy Surveys DR9 Tractor forced WISE photometry (cross-check) 5. radio CDS XMatch: FIRST (VIII/92/first14) and NVSS (VIII/65/nvss) 5b. galex CDS XMatch: GALEX GR6+7 AIS (II/335/galex_ais), NUV/FUV 6. wise_rsr WISE relative spectral response curves (IPAC) 7. brown14 Brown et al. 2014 SED atlas (CDS J/ApJS/212/18): table1 + spectra of the early-type galaxies """ import io import os import subprocess import sys import time import urllib.parse import numpy as np import pandas as pd import mirq_common as M SKY = "https://skyserver.sdss.org/dr17/SkyServerWS/SearchTools/SqlSearch" IRSA = "https://irsa.ipac.caltech.edu/TAP/sync" DATALAB = "https://datalab.noirlab.edu/tap/sync" XMATCH = "http://cdsxmatch.u-strasbg.fr/xmatch/api/v1/sync" CDS_B14 = "https://cdsarc.cds.unistra.fr/ftp/J/ApJS/212/18/" RSR_URL = "https://wise2.ipac.caltech.edu/docs/release/prelim/expsup/figures/RSR-W%d.txt" SAMPLE = os.path.join(M.DATA, "sample_etg.csv.gz") SAMPLE_COLS = ["id", "specObjID", "plate", "mjd", "fiber", "ra", "dec", "z", "zhelio", "veldisp", "e_logs", "theta_e", "logr", "logi", "Mr", "msel", "devmag_r", "snr", "conc", "ba", "fracdev", "gr", "petroR50", "devrad", "ew_ha_port", "ew_ha_mpa", "d4000", "bpt", "clean", "nchild", "ext_r"] FP_CACHE = "fermi_raw/fp_kiii/derived/sample_sdss.csv" # ---------------------------------------------------------------------------------------------------- def get_sample(): if os.path.exists(SAMPLE): return pd.read_csv(SAMPLE) src = FP_CACHE if not os.path.exists(src): # rebuild with the fp_kiii code (SDSS SkyServer, ~1 h) env = dict(os.environ, FP_KIII_RAW=os.path.join(M.RAW, "fp_kiii_rebuild")) bdir = os.path.join(M.HERE, "fp_kiii_build") for script in ("download_data.py", "build_samples.py"): subprocess.check_call([sys.executable, os.path.join(bdir, script)], env=env, cwd=bdir) src = os.path.join(env["FP_KIII_RAW"], "derived", "sample_sdss.csv") d = pd.read_csv(src) d = d[d.in_sample == 1][SAMPLE_COLS] d.to_csv(SAMPLE, index=False, float_format="%.7g", compression="gzip") print("sample written: %d rows" % len(d)) return pd.read_csv(SAMPLE) # ---------------------------------------------------------------------------------------------------- LINE_SQL = """ SELECT s.specObjID, x.bptclass, x.lgm_tot_p50, x.sfr_tot_p50, x.specsfr_tot_p50, l.h_alpha_flux, l.h_alpha_flux_err, l.h_beta_flux, l.h_beta_flux_err, l.oiii_5007_flux, l.oiii_5007_flux_err, l.nii_6584_flux, l.nii_6584_flux_err, l.h_alpha_eqw, l.h_alpha_eqw_err, l.oii_3729_eqw, i.d4000_n, i.lick_hd_a, i.lick_hd_a_err, p.modelMag_u, p.modelMag_g, p.modelMag_r, p.modelMag_i, p.modelMag_z, p.cModelMag_r, p.extinction_u, p.extinction_g, p.extinction_r, p.extinction_i, p.extinction_z FROM SpecObj s JOIN PhotoObjAll p ON s.bestObjID = p.objID LEFT JOIN galSpecExtra x ON x.specObjID = s.specObjID LEFT JOIN galSpecLine l ON l.specObjID = s.specObjID LEFT JOIN galSpecIndx i ON i.specObjID = s.specObjID WHERE s.class='GALAXY' AND s.zWarning=0 AND s.sciencePrimary=1 AND s.survey='sdss' AND s.z BETWEEN 0.0033 AND 0.125 AND s.velDisp > 50 AND p.fracDeV_r >= 0.5 AND s.ra >= {lo} AND s.ra < {hi} """ def get_sdss_lines(step=4): dest_all = os.path.join(M.RAW, "sdss_lines.csv") if os.path.exists(dest_all): return parts = [] for lo in range(0, 360, step): dest = os.path.join(M.RAW, "sdss_lines", "lines_ra%03d.csv" % lo) sql = " ".join(LINE_SQL.format(lo=lo, hi=lo + step).split()) url = SKY + "?" + urllib.parse.urlencode({"cmd": sql, "format": "csv"}) def chk(blob): if b"specObjID" not in blob[:200]: raise RuntimeError("bad payload %r" % blob[:200]) M.fetch(url, dest, timeout=600, min_bytes=100, check=chk) x = pd.read_csv(dest, comment="#") if len(x): parts.append(x) print(" sdss lines chunk %3d: %d rows" % (lo, len(x)), flush=True) pd.concat(parts, ignore_index=True).drop_duplicates("specObjID").to_csv(dest_all, index=False) # ---------------------------------------------------------------------------------------------------- AW_COLS = ("designation ra dec w1mpro w1sigmpro w1snr w1rchi2 w2mpro w2sigmpro w2snr w2rchi2 w3mpro w3sigmpro w3snr " "w3rchi2 w4mpro w4sigmpro w4snr w4rchi2 rchi2 nb na w1sat w3sat w4sat cc_flags ext_flg ph_qual det_bit " "moon_lev w1cov w2cov w3cov w4cov w1flux w1sigflux w2flux w2sigflux w3flux w3sigflux w4flux w4sigflux " "w3sky w3sigsk w3conf w4sky w4sigsk w4conf " "w1mag w1sigm w1flg w2mag w2sigm w2flg w3mag w3sigm w3flg w4mag w4sigm w4flg " "w1mag_4 w1sigm_4 w1flg_4 w3mag_4 w3sigm_4 w3flg_4 w4mag_4 w4sigm_4 w4flg_4 " "w1mag_8 w1sigm_8 w1flg_8 w3mag_8 w3sigm_8 w3flg_8 w4mag_8 w4sigm_8 w4flg_8 " "xscprox w1rsemi w1ba w1gmag w1gerr w1gflg w2gmag w2gerr w2gflg w3rsemi w3gmag w3gerr w3gflg " "w4rsemi w4gmag w4gerr w4gflg tmass_key r_2mass n_2mass k_m_2mass").split() def get_allwise(samp, chunk=5000, radius_arcsec=6.0): dest_all = os.path.join(M.RAW, "allwise_match.csv") if os.path.exists(dest_all): return cols = ",".join("a." + c for c in AW_COLS) q = ("SELECT u.sid,%s FROM allwise_p3as_psd a, TAP_UPLOAD.up u WHERE " "CONTAINS(POINT('ICRS',a.ra,a.dec),CIRCLE('ICRS',u.ra_in,u.dec_in,%.6f))=1" % (cols, radius_arcsec / 3600.0)) parts = [] for i0 in range(0, len(samp), chunk): dest = os.path.join(M.RAW, "allwise", "aw_%06d.csv" % i0) sub = samp.iloc[i0:i0 + chunk] up = pd.DataFrame({"sid": sub.id.values, "ra_in": sub.ra.values, "dec_in": sub["dec"].values}) blob = up.to_csv(index=False, float_format="%.7f").encode() body, hdr = M.multipart(dict(UPLOAD="up,param:table.csv", FORMAT="CSV", QUERY=q), {"table.csv": ("table.csv", blob)}) def chk(b): if not b.startswith(b"sid,designation"): raise RuntimeError("bad payload %r" % b[:300]) M.fetch(IRSA, dest, data=body, headers=hdr, timeout=1200, check=chk) x = pd.read_csv(dest, dtype={"cc_flags": str, "ph_qual": str}) parts.append(x) print(" allwise chunk %6d: %d matches for %d positions" % (i0, len(x), len(sub)), flush=True) pd.concat(parts, ignore_index=True).to_csv(dest_all, index=False) # ---------------------------------------------------------------------------------------------------- LS_Q = ("SELECT s.specobjid, x.distance, t.ls_id, t.type, t.shape_r, t.sersic, t.flux_g, t.flux_r, t.flux_z, " "t.flux_w1, t.flux_w2, t.flux_w3, t.flux_w4, t.flux_ivar_w1, t.flux_ivar_w2, t.flux_ivar_w3, t.flux_ivar_w4, " "t.fracflux_w1, t.fracflux_w3, t.fracflux_w4, t.rchisq_w1, t.rchisq_w3, t.rchisq_w4, t.nobs_w3, t.nobs_w4, " "t.mw_transmission_r, t.maskbits, t.fracflux_r " "FROM sdss_dr17.specobj s JOIN sdss_dr17.x1p5__specobj__ls_dr9__tractor x ON x.id1=s.specobjid " "JOIN ls_dr9.tractor t ON t.ls_id=x.id2 " "WHERE s.ra >= {lo} AND s.ra < {hi} AND s.z BETWEEN 0.015 AND 0.125 AND s.class='GALAXY' AND s.veldisp > 50") def get_lsdr9(samp, step=4): dest_all = os.path.join(M.RAW, "lsdr9_match.csv") if os.path.exists(dest_all): return ids = set(samp.specObjID.values) parts = [] for lo in range(0, 360, step): dest = os.path.join(M.RAW, "lsdr9", "ls_ra%03d.csv" % lo) data = urllib.parse.urlencode(dict(REQUEST="doQuery", LANG="ADQL", FORMAT="csv", QUERY=LS_Q.format(lo=lo, hi=lo + step))).encode() def chk(b): if not b.startswith(b"specobjid,distance"): raise RuntimeError("bad payload %r" % b[:300]) M.fetch(DATALAB, dest, data=data, timeout=1200, min_bytes=20, check=chk) x = pd.read_csv(dest) x = x[x.specobjid.isin(ids)] parts.append(x) print(" ls dr9 chunk %3d: %d sample rows" % (lo, len(x)), flush=True) out = pd.concat(parts, ignore_index=True).sort_values("distance").drop_duplicates("specobjid") out.to_csv(dest_all, index=False) # ---------------------------------------------------------------------------------------------------- def get_radio(samp): for cat, name, rad in (("vizier:VIII/92/first14", "first", 10.0), ("vizier:VIII/65/nvss", "nvss", 30.0)): dest = os.path.join(M.RAW, "xmatch_%s.csv" % name) if os.path.exists(dest): continue parts = [] for i0 in range(0, len(samp), 40000): part = os.path.join(M.RAW, "radio", "%s_%06d.csv" % (name, i0)) blob = samp.iloc[i0:i0 + 40000][["id", "ra", "dec"]].to_csv(index=False).encode() body, hdr = M.multipart(dict(request="xmatch", distMaxArcsec=str(rad), RESPONSEFORMAT="csv", colRA1="ra", colDec1="dec", cat2=cat), {"cat1": ("cat1.csv", blob)}) def chk(b): if not b.startswith(b"angDist"): raise RuntimeError("bad payload %r" % b[:300]) M.fetch(XMATCH, part, data=body, headers=hdr, timeout=1200, min_bytes=20, check=chk) parts.append(pd.read_csv(part)) out = pd.concat(parts, ignore_index=True) out.to_csv(dest, index=False) print(" %s: %d matches within %.0f arcsec" % (name, len(out), rad), flush=True) def get_galex(samp): """GALEX GR6+7 AIS (Bianchi+ 2017, VizieR II/335/galex_ais), 4 arcsec: NUV for the residual-star-formation test.""" dest = os.path.join(M.RAW, "xmatch_galex.csv") if os.path.exists(dest): return parts = [] for i0 in range(0, len(samp), 20000): part = os.path.join(M.RAW, "galex", "galex_%06d.csv" % i0) blob = samp.iloc[i0:i0 + 20000][["id", "ra", "dec"]].to_csv(index=False).encode() body, hdr = M.multipart(dict(request="xmatch", distMaxArcsec="4", RESPONSEFORMAT="csv", colRA1="ra", colDec1="dec", cat2="vizier:II/335/galex_ais"), {"cat1": ("cat1.csv", blob)}) def chk(b): if not b.startswith(b"angDist"): raise RuntimeError("bad payload %r" % b[:300]) M.fetch(XMATCH, part, data=body, headers=hdr, timeout=1200, min_bytes=20, check=chk) x = pd.read_csv(part, usecols=["angDist", "id", "FUVmag", "e_FUVmag", "NUVmag", "e_NUVmag", "E(B-V)", "NUVexp", "FUVexp"]) parts.append(x) out = pd.concat(parts, ignore_index=True) out.to_csv(dest, index=False) print(" galex: %d matches within 4 arcsec" % len(out), flush=True) # ---------------------------------------------------------------------------------------------------- def get_templates(): for b in (1, 2, 3, 4): M.fetch(RSR_URL % b, os.path.join(M.RAW, "wise_rsr", "RSR-W%d.txt" % b), min_bytes=1000) t1 = M.fetch(CDS_B14 + "table1.dat", os.path.join(M.RAW, "brown14", "table1.dat"), min_bytes=1000) M.fetch(CDS_B14 + "table3.dat", os.path.join(M.RAW, "brown14", "table3.dat"), min_bytes=1000) M.fetch(CDS_B14 + "ReadMe", os.path.join(M.RAW, "brown14", "ReadMe"), min_bytes=1000) n = 0 for line in open(t1): name = line[:15].strip() try: ttype = float(line[49:53]) except ValueError: continue if ttype <= -1.0: # E and S0 galaxies of the atlas fn = name.replace(" ", "_") + "_spec.dat" M.fetch(CDS_B14 + "sp/" + urllib.parse.quote(fn), os.path.join(M.RAW, "brown14", "sp", fn), min_bytes=1000) n += 1 print(" Brown+14 early-type spectra: %d" % n, flush=True) def main(): t0 = time.time() samp = get_sample() print("sample: %d galaxies" % len(samp), flush=True) get_templates() get_radio(samp) get_galex(samp) get_allwise(samp) get_sdss_lines() get_lsdr9(samp) print("downloads complete (%.0f s)" % (time.time() - t0), flush=True) if __name__ == "__main__": main()