#!/usr/bin/env python """Build a cross-band merged Vega-mag catalog for wd1 / wd2 (EWOCS comparison). wd1/wd2 have only per-FILTER merged catalogs (merged across exposures), not a cross-FILTER merged multiband table like w51 already has. This does a positional UNION cross-match of the per-filter `_vetted` products and attaches per-filter Vega magnitudes computed with the SAME recipe the jwst-gc-pipeline uses: flux_jy = flux[MJy/sr] * pixel_area[sr] (pixel scale from cat meta PIXSCALE) mag_vega = -2.5 * log10(flux_jy / ZP_svo_Jy) (ZP from SVO Filter Profile Service) NOTE (scope): this is a UNION of per-filter detections matched within MATCH_RADIUS. It is NOT the full pipeline merge -- no cross-program forced photometry, no satstar wing-cal reconciliation. Adequate for comparison CMD/CCDs; documented as such. """ import os import re import glob import warnings import numpy as np import astropy.units as u from astropy.table import Table, Column from astropy.coordinates import SkyCoord from astroquery.svo_fps import SvoFps MATCH_RADIUS = 0.10 * u.arcsec BASE = "/orange/adamginsburg/jwst" OUTDIR = f"{BASE}/w51-ewocs/catalogs" # reference filter for the union seed: a deep broadband present in both fields. REF_FILTER = "f200w" _JF = None def svo_zp(filtername): """Vega zeropoint [Jy] for a NIRCam filter from the SVO service.""" global _JF if _JF is None: _JF = SvoFps.get_filter_list("JWST") _JF.add_index("filterID") fid = f"JWST/NIRCam.{filtername.upper()}" return u.Quantity(_JF.loc[fid]["ZeroPoint"], u.Jy) def best_catalog(field, filt): """Deepest merged (not per-module) resbgsub vetted catalog for a filter.""" cands = glob.glob(f"{BASE}/{field}/catalogs/{filt}_merged_indivexp_merged_*" f"_dao_basic_vetted.fits") cands = [c for c in cands if "nrca" not in c and "nrcb" not in c] if not cands: return None def rank(path): b = os.path.basename(path) m = re.search(r"_m(\d+)_", b) return ("resbgsub" in b, int(m.group(1)) if m else -1) return sorted(cands, key=rank)[-1] def load_filter(field, filt): """Load one filter's vetted catalog, return (SkyCoord, per-filter columns dict).""" path = best_catalog(field, filt) if path is None: return None t = Table.read(path) px = t.meta["PIXSCALE"] # arcsec/pixel pixarea = ((px * u.arcsec) ** 2).to(u.sr) flux = np.asarray(t["flux"], float) ferr = np.asarray(t["flux_err_prop" if "flux_err_prop" in t.colnames else "flux_err"], float) with np.errstate(all="ignore"): flux_jy = (flux * u.MJy / u.sr * pixarea).to(u.Jy) vega = (-2.5 * np.log10(flux_jy / svo_zp(filt))).value emag = 2.5 / np.log(10) * np.abs(ferr / flux) sc = SkyCoord(t["skycoord"]) cols = { f"mag_vega_{filt}": vega, f"emag_ab_{filt}": emag, f"qfit_{filt}": np.asarray(t["qfit"], float), f"is_saturated_{filt}": np.asarray(t["is_saturated"], bool), f"flux_jy_{filt}": np.asarray(flux_jy.value, float), } print(f" {filt}: {len(t)} src (PIXSCALE={px:.5f}\", ZP={svo_zp(filt):.1f}, " f"vega med={np.nanmedian(vega):.2f}) <- {os.path.basename(path)}") return sc, cols def build(field, filters): print(f"=== building merged catalog for {field} ===") loaded = {} for f in filters: r = load_filter(field, f) if r is not None: loaded[f] = r present = [f for f in filters if f in loaded] if REF_FILTER in present: order = [REF_FILTER] + [f for f in present if f != REF_FILTER] else: order = present print(f" union order (ref first): {order}") # --- incremental UNION build --- # master_sc: growing SkyCoord of union positions (seeded by ref filter). sc0, cols0 = loaded[order[0]] n = len(sc0) master_ra = list(sc0.ra.deg) master_dec = list(sc0.dec.deg) # per-filter column store, indexed by master row; fill masked defaults store = {} def blank(colname, template): arr = np.full(len(master_ra), np.nan) if template.dtype == bool: arr = np.zeros(len(master_ra), bool) return arr # seed columns from ref filter for k, v in cols0.items(): arr = np.array(v, dtype=v.dtype) store[k] = arr for f in order[1:]: sc, cols = loaded[f] master_sc = SkyCoord(master_ra * u.deg, master_dec * u.deg) idx, sep, _ = sc.match_to_catalog_sky(master_sc) matched = sep < MATCH_RADIUS # ensure store has columns for this filter (all-blank then fill) for k, v in cols.items(): if k not in store: store[k] = (np.zeros(len(master_ra), bool) if v.dtype == bool else np.full(len(master_ra), np.nan)) # fill matched for k, v in cols.items(): store[k][idx[matched]] = v[matched] # append unmatched as new union rows newsrc = ~matched n_new = int(newsrc.sum()) if n_new: master_ra.extend(sc.ra.deg[newsrc]) master_dec.extend(sc.dec.deg[newsrc]) for k in store: pad = (np.zeros(n_new, bool) if store[k].dtype == bool else np.full(n_new, np.nan)) store[k] = np.concatenate([store[k], pad]) for k, v in cols.items(): store[k][-n_new:] = v[newsrc] print(f" +{f}: matched {int(matched.sum())}, new {n_new}, " f"union now {len(master_ra)}") out = Table() out["skycoord_ref"] = SkyCoord(master_ra * u.deg, master_dec * u.deg) for k in sorted(store): out[k] = store[k] out.meta["FIELD"] = field out.meta["FILTERS"] = ",".join(order) out.meta["NOTE"] = ("UNION positional cross-match of per-filter vetted DAO " "catalogs; Vega mags via SVO ZP + PIXSCALE. Not the full " "pipeline merge.") outfn = f"{OUTDIR}/{field}_ewocs_merged_vega.fits" out.write(outfn, overwrite=True) print(f" wrote {outfn} ({len(out)} union sources, {len(out.colnames)} cols)") return out WD1 = ["f115w", "f150w", "f164n", "f187n", "f200w", "f212n", "f277w", "f323n", "f405n", "f444w", "f466n"] WD2 = ["f115w", "f150w", "f162m", "f164n", "f182m", "f187n", "f200w", "f212n", "f250m", "f277w", "f300m", "f323n", "f335m", "f405n", "f410m", "f444w", "f466n"] if __name__ == "__main__": with warnings.catch_warnings(): warnings.simplefilter("ignore") build("wd1", WD1) build("wd2", WD2)