"""Step 9: independent aperture photometry of candidates and controls on Legacy Survey r,z images. Why: the ALFALFA-SDSS catalogue uses SDSS pipeline cmodel magnitudes (and unWISE PSF-based fluxes), which are known to shred or truncate large/low-surface-brightness galaxies. A galaxy whose catalogue flux is too faint by >~1 mag looks exactly like a "harvested" galaxy on the Tully-Fisher diagram. Method (deliberately simple, no external packages): 1. sky = median of the outer frame (r > 120 px), per band; noise from MAD. 2. detection image = r+z smoothed with a 1.5 px Gaussian; segmentation at 2.5 sigma (scipy.ndimage.label). 3. target = segment under the image centre (or nearest within 8 px). 4. flux-weighted second moments of the target segment -> centre, axis ratio, PA, rms semi-major axis. 5. elliptical curve of growth in 2 px annuli; pixels of OTHER segments are masked and replaced by the annulus median. Two fluxes: 'mean' (annulus mean x area; includes superposed stars) and 'med' (annulus median x area; robust against superposed stars, slightly low for clumpy discs). The aperture is 4 x the rms semi-major axis (>= 8"), capped at 140". 6. zero points are calibrated empirically on control galaxies that have Cosmicflows-4 total magnitudes. 7. (revision) the measurement is only USED when the photometered segment is the target itself: ap_valid = segment centroid within CEN_MAX pixels of the catalogue position and at most FM_MAX of the aperture masked. Without this check the "target" can be the halo of a nearby bright star or a neighbouring galaxy (referee point E1: AGC 320380, 112351, 103316). Output: data/aperture_phot.csv """ import os import warnings import numpy as np import pandas as pd from scipy import ndimage from astropy.io import fits from common import RAW, DATA, mad_sigma warnings.filterwarnings("ignore") CEN_MAX = 6.0 # pixels = arcsec FM_MAX = 0.20 def photometer(path): with fits.open(path) as h: cube = np.asarray(h[0].data, float) if cube.ndim != 3 or cube.shape[0] < 2: return None ny, nx = cube.shape[1:] yy, xx = np.mgrid[0:ny, 0:nx] cy, cx = (ny - 1) / 2.0, (nx - 1) / 2.0 rr = np.hypot(yy - cy, xx - cx) outer = rr > 120 out = {} imgs = [] for b, name in enumerate("rz"): im = cube[b] if not np.isfinite(im).any() or np.all(im == 0): return None sky = np.nanmedian(im[outer]) for _ in range(3): # clip bright sources out of the sky estimate sel = outer & (np.abs(im - sky) < 3 * mad_sigma(im[outer & np.isfinite(im)])) sky = np.nanmedian(im[sel]) imgs.append(np.nan_to_num(im - sky)) out["sky_sig_" + name] = mad_sigma(im[outer & np.isfinite(im)]) det = ndimage.gaussian_filter(imgs[0] + imgs[1], 1.5) sig = mad_sigma(det[outer]) if not sig > 0: return None seg, nlab = ndimage.label(det > 2.5 * sig) icy, icx = int(round(cy)), int(round(cx)) lab = seg[icy, icx] if lab == 0: sub = seg[icy - 8:icy + 9, icx - 8:icx + 9] d = np.hypot(*(np.mgrid[-8:9, -8:9])) cand = [(d[sub == l].min(), l) for l in np.unique(sub) if l > 0] lab = min(cand)[1] if cand else 0 out["detected"] = bool(lab > 0) if lab > 0: m = seg == lab w = np.clip(det[m], 0, None) x0 = np.sum(xx[m] * w) / w.sum() y0 = np.sum(yy[m] * w) / w.sum() dx, dy = xx[m] - x0, yy[m] - y0 sxx, syy, sxy = np.sum(w * dx * dx) / w.sum(), np.sum(w * dy * dy) / w.sum(), np.sum(w * dx * dy) / w.sum() tr, dt = sxx + syy, np.sqrt(max(((sxx - syy) / 2) ** 2 + sxy ** 2, 0)) a2, b2 = tr / 2 + dt, max(tr / 2 - dt, 0.5) a_rms, q = np.sqrt(a2), np.sqrt(b2 / a2) pa = 0.5 * np.arctan2(2 * sxy, sxx - syy) npix = int(m.sum()) else: x0, y0, a_rms, q, pa, npix = cx, cy, 3.0, 1.0, 0.0, 0 q = float(np.clip(q, 0.15, 1.0)) r_ap = float(np.clip(4.0 * a_rms, 8.0, 140.0)) ct, st = np.cos(pa), np.sin(pa) u = (xx - x0) * ct + (yy - y0) * st v = -(xx - x0) * st + (yy - y0) * ct rell = np.hypot(u, v / q) other = ndimage.binary_dilation((seg > 0) & (seg != lab), iterations=2) edges = np.arange(0, r_ap + 2.0, 2.0) for b, name in enumerate("rz"): im = imgs[b] fmean = fmed = 0.0 for lo, hi in zip(edges[:-1], edges[1:]): ann = (rell >= lo) & (rell < hi) good = ann & ~other if good.sum() < 3: continue fmean += np.mean(im[good]) * ann.sum() fmed += np.median(im[good]) * ann.sum() out["fmean_" + name], out["fmed_" + name] = fmean, fmed out.update(x0=x0, y0=y0, a_rms=a_rms, q_img=q, pa=np.degrees(pa), r_ap=r_ap, seg_npix=npix, cen_off=float(np.hypot(x0 - cx, y0 - cy)), frac_masked=float(other[rell < r_ap].mean())) return out def main(): t = pd.read_csv(os.path.join(DATA, "vet_targets.csv")) rows = [] for r in t.itertuples(): p = os.path.join(RAW, "cutouts", f"ls_{r.uid}.fits") if not (r.want_fits and os.path.exists(p) and os.path.getsize(p) > 500): continue try: o = photometer(p) except Exception as e: # corrupt download etc. o = None if o is None: rows.append(dict(uid=r.uid, phot_ok=False)) continue o.update(uid=r.uid, phot_ok=True) rows.append(o) a = pd.DataFrame(rows) a = a.merge(t[["uid", "src", "role", "PGC", "AGC", "D", "m_i", "V"]], on="uid", how="left") # ---- empirical zero points from control galaxies with CF4 total (raw) r and z magnitudes cf = pd.read_csv(os.path.join(RAW, "cf4_t1.csv"), usecols=["PGC", "rmag", "imag", "zmag", "QSflag"]) cf = cf.rename(columns={"rmag": "cf_r", "imag": "cf_i", "zmag": "cf_z"}) a = a.merge(cf, on="PGC", how="left") for name in "rz": for kind in ("mean", "med"): f = a[f"f{kind}_{name}"] a[f"minst_{kind}_{name}"] = -2.5 * np.log10(f.where(f > 0)) cal = a[(a.role == "control") & (a.QSflag >= 4) & a.phot_ok & a.detected.fillna(False).astype(bool) & a.cf_r.between(8, 19) & a.cf_z.between(8, 19)] zp = {} for name in "rz": for kind in ("mean", "med"): d = cal[f"cf_{name}"] - cal[f"minst_{kind}_{name}"] zp[(kind, name)] = (d.median(), mad_sigma(d.dropna()), int(d.notna().sum())) a[f"m_{kind}_{name}"] = a[f"minst_{kind}_{name}"] + d.median() print(" zero points vs CF4 total magnitudes (median, robust scatter, N):") for k, v in zp.items(): print(f" {k[0]:4s} {k[1]}: ZP={v[0]:.3f} scatter={v[1]:.3f} mag N={v[2]}") # ---- synthetic i-band: mean of r and z plus a colour offset calibrated on CF4 i (raw total) for kind in ("mean", "med"): rz = 0.5 * (a[f"m_{kind}_r"] + a[f"m_{kind}_z"]) c = (a.cf_i - rz)[cal.index] a[f"i_ap_{kind}"] = rz + c.median() print(f" synthetic i ({kind}): offset {c.median():+.3f}, scatter vs CF4 i {mad_sigma(c.dropna()):.3f} mag") # ---- is the photometered segment the target? a["ap_valid"] = (a.phot_ok & a.detected.fillna(False).astype(bool) & (a.cen_off <= CEN_MAX) & (a.frac_masked <= FM_MAX)) for role in ("candidate", "control", "control_ap"): m = (a.role == role) & a.phot_ok print(f" segment centred on the target (<= {CEN_MAX:.0f} px, masked <= {FM_MAX}): {role:10s} " f"{a.ap_valid[m].sum()}/{m.sum()}") # ---- catalogue-minus-aperture magnitude for ALFALFA-only galaxies (SDSS cmodel i, uncorrected) al = a.src == "alfalfa" ctrl = al & a.role.isin(["control", "control_ap"]) & a.ap_valid for kind in ("mean", "med"): d = a.m_i - a[f"i_ap_{kind}"] off = d[ctrl].median() a[f"dmag_{kind}"] = np.where(al, d - off, np.nan) # >0: catalogue fainter than the image flux print(f" ALFALFA-only controls: cmodel-aperture ({kind}) median {off:+.3f}, robust scatter " f"{mad_sigma(d[ctrl].dropna()):.3f} mag, N={ctrl.sum()}; fraction with catalogue >0.75 mag too faint: " f"{(d[ctrl] - off > 0.75).mean():.3f}") a.to_csv(os.path.join(DATA, "aperture_phot.csv"), index=False) print(f" aperture photometry: {a.phot_ok.sum()} measured, {(~a.phot_ok).sum()} without usable image") if __name__ == "__main__": main()