"""Step 8: download per-object vetting data for union-sample candidates and a random control sample. Targets : (a) every union galaxy with D < log10(0.5) (the f>=0.75 and f>=0.9 candidates are subsets), (b) a seeded random control sample of 300 union non-candidates, (c) 200 further seeded random ALFALFA-only non-candidates ("control_ap"; FITS cutout only), used solely to measure how often the aperture-photometry check of s09 would wrongly excuse a galaxy. Per target: * SDSS DR16 cone (VizieR V/154, r = 2.5 arcmin): redshifts of other galaxies inside the Arecibo beam * Legacy Survey r,z FITS cutout (5 arcmin, 1"/pix) for independent aperture photometry (ALFALFA-only galaxies, whose catalogue photometry is SDSS-pipeline cmodel, and all controls) * JPG cutouts (wide 4.3 arcmin and zoom 1.3 arcmin) for the f>=0.75 candidates and 40 controls Output : data/vet_targets.csv ; files under RAW/cones and RAW/cutouts """ import os import sys import numpy as np import pandas as pd from concurrent.futures import ThreadPoolExecutor from common import RAW, DATA, VIZIER_TAP, tap_csv, http_get T50, T75, T90 = np.log10(0.5), np.log10(0.25), np.log10(0.1) N_CONTROL = 300 N_CONTROL_VISUAL = 40 N_CONTROL_AP = 200 SEED = 20261011 SEED_AP = 20261012 CONE_DEG = 2.5 / 60.0 VIEWER = "https://www.legacysurvey.org/viewer/" def targets(): u = pd.read_csv(os.path.join(DATA, "union_sample.csv")) cand = u[u.D < T50].copy() cand["role"] = "candidate" cand["ctrl_order"] = -1 rng = np.random.default_rng(SEED) pool = u[u.D >= T50] ctrl = pool.iloc[rng.permutation(len(pool))[:N_CONTROL]].copy() ctrl["role"] = "control" ctrl["ctrl_order"] = np.arange(len(ctrl)) # extra aperture-photometry controls: appended AFTER the original targets so that their order is unchanged rng2 = np.random.default_rng(SEED_AP) pool2 = u[(u.D >= T50) & (u.src == "alfalfa") & ~u.uid.isin(ctrl.uid)] cap = pool2.iloc[rng2.permutation(len(pool2))[:N_CONTROL_AP]].copy() cap["role"] = "control_ap" cap["ctrl_order"] = np.arange(len(cap)) t = pd.concat([cand, ctrl, cap], ignore_index=True) t["visual"] = ((t.role == "candidate") & (t.D < T75)) | ((t.role == "control") & (t.ctrl_order < N_CONTROL_VISUAL)) t["want_fits"] = (t.src == "alfalfa") | (t.role == "control") | t.visual t["want_cone"] = t.role != "control_ap" return t def fetch_cone(uid, ra, dec): q = ('SELECT objID, RA_ICRS, DE_ICRS, "class", "mode", gmag, rmag, imag, zsp, zph, e_zph ' 'FROM "V/154/sdss16" WHERE 1=CONTAINS(POINT(\'ICRS\',RA_ICRS,DE_ICRS),' f'CIRCLE(\'ICRS\',{ra:.6f},{dec:.6f},{CONE_DEG:.5f})) AND "class"=3 AND rmag<20.5 AND "mode"=1') path = os.path.join(RAW, f"cones/sdss16_{uid}.csv") for attempt in range(4): tap_csv(VIZIER_TAP, q, f"cones/sdss16_{uid}.csv", force=attempt > 0) with open(path) as f: if f.readline().startswith("objID"): # header present (an empty cone is header-only) return path raise RuntimeError("cone query returned no header") def fetch_file(url, path, minsize=500, retries=3, timeout=120): if os.path.exists(path) and os.path.getsize(path) > minsize: return path blob = http_get(url, timeout=timeout, retries=retries) with open(path, "wb") as f: f.write(blob) return path def fetch_fits(uid, ra, dec): """r,z cutout; DR9 first (covers the ALFALFA sky), then DR10 (southern sky).""" path = os.path.join(RAW, "cutouts", f"ls_{uid}.fits") if os.path.exists(path) and os.path.getsize(path) > 500: return path if os.path.exists(path + ".missing"): # both layers failed in an earlier run (no imaging there) return None last = None for layer in ("ls-dr9", "ls-dr10"): try: # where a layer has no imaging the server answers HTTP 500, so do not retry many times return fetch_file(f"{VIEWER}cutout.fits?ra={ra:.5f}&dec={dec:.5f}&layer={layer}&pixscale=1.0&size=300&bands=rz", path, retries=2, timeout=60) except Exception as e: # noqa last = e open(path + ".missing", "w").write(str(last)[:200]) raise RuntimeError(last) def ps1_jpg(ra, dec, size_arcsec, path): """Pan-STARRS1 g,r,i colour cutout (fallback where the Legacy Surveys and SDSS have no imaging).""" tab = http_get(f"https://ps1images.stsci.edu/cgi-bin/ps1filenames.py?ra={ra:.5f}&dec={dec:.5f}&filters=gri", timeout=60).decode().strip().splitlines() names = {ln.split()[4]: ln.split()[7] for ln in tab[1:]} if not all(b in names for b in "gri"): raise RuntimeError("no PS1 stack") npix = int(round(size_arcsec / 0.25)) url = (f"https://ps1images.stsci.edu/cgi-bin/fitscut.cgi?red={names['i']}&green={names['r']}&blue={names['g']}" f"&ra={ra:.5f}&dec={dec:.5f}&size={npix}&output_size=256&format=jpg&autoscale=99.5") return fetch_file(url, path) def fetch_jpgs(uid, ra, dec): for tag, pix in (("wide", 1.0), ("zoom", 0.3)): path = os.path.join(RAW, "cutouts", f"{tag}_{uid}.jpg") done = False for layer in ("ls-dr9", "ls-dr10", "sdss"): lpath = path if layer == "ls-dr9" else path.replace(".jpg", f"_{layer}.jpg") try: fetch_file(f"{VIEWER}cutout.jpg?ra={ra:.5f}&dec={dec:.5f}&layer={layer}&pixscale={pix}&size=256", lpath) except Exception: continue if os.path.getsize(lpath) > 6000: # an empty (no-coverage) JPG is only a few kB done = True break if not done: ps1_jpg(ra, dec, 256 * pix, path.replace(".jpg", "_ps1.jpg")) def main(): t = targets() t.to_csv(os.path.join(DATA, "vet_targets.csv"), index=False) print(f" targets: {len(t)} ({(t.role == 'candidate').sum()} candidates, {(t.role == 'control').sum()} controls, " f"{(t.role == 'control_ap').sum()} extra aperture-photometry controls;" f" {t.visual.sum()} for visual inspection; {t.want_fits.sum()} FITS cutouts)") def work(r): nbad = 0 for fn, args, want in ((fetch_cone, (r.uid, r.ra, r.dec), r.want_cone), (fetch_fits, (r.uid, r.ra, r.dec), r.want_fits), (fetch_jpgs, (r.uid, r.ra, r.dec), r.visual)): if not want: continue try: fn(*args) except Exception as e: # keep going; missing files are handled downstream nbad += 1 print(f" WARN {r.uid} {fn.__name__}: {str(e)[:150]}", file=sys.stderr) return nbad # network-bound, so a few threads are used (no heavy CPU work here) with ThreadPoolExecutor(max_workers=4) as ex: nfail = sum(ex.map(work, list(t.itertuples()))) print(f" vetting downloads finished, failed items: {nfail}") if __name__ == "__main__": main()