diff --git a/autoreduce/instruments/__init__.py b/autoreduce/instruments/__init__.py index 0704c51..2cd02f7 100644 --- a/autoreduce/instruments/__init__.py +++ b/autoreduce/instruments/__init__.py @@ -7,3 +7,5 @@ from .adapter import InstrumentAdapter, get, register, registered_keys from .acs_wfc import ACS_WFC +from .wfc3_uvis import WFC3_UVIS +from .wfc3_ir import WFC3_IR diff --git a/autoreduce/instruments/adapter.py b/autoreduce/instruments/adapter.py index 2f33d9c..8bcd37b 100644 --- a/autoreduce/instruments/adapter.py +++ b/autoreduce/instruments/adapter.py @@ -20,7 +20,10 @@ class InstrumentAdapter: crds_reference_subpath: str # where CRDS syncs this instrument's files supports_cte_correction: bool default_drizzle_kwargs: Dict[str, object] - saturation_dn: float # conservative full-well / saturation level + saturation_dn: float # conservative full-well / saturation level, electrons + # The adapter's recommendation for TargetSpec.final_scale (which remains + # the user-facing dial); documents sensible sampling for this detector. + recommended_final_scale: float = 0.05 def scale_ratio(self, final_scale: float) -> float: """s = output scale / native scale, as used by the Casertano factor.""" diff --git a/autoreduce/instruments/wfc3_ir.py b/autoreduce/instruments/wfc3_ir.py new file mode 100644 index 0000000..5ea229a --- /dev/null +++ b/autoreduce/instruments/wfc3_ir.py @@ -0,0 +1,34 @@ +""" +WFC3/IR — adapter #3 (roadmap phase 2), the genuinely different path. + +No CTE correction exists for the IR channel, so the calibrated product is +``_flt`` (already in e-/s from up-the-ramp fitting, which also rejects most +cosmic rays per read — AstroDrizzle's driz_cr then handles the residue when +multiple exposures exist). Native scale ~0.128″/pix under-samples the PSF, so +dithered programs conventionally drizzle to a finer grid; 0.065″/pix is the +adapter's recommendation (half-native, within the 0.06–0.08 range common in +deep-field practice) — `TargetSpec.final_scale` remains the user dial and +star-poor or poorly-dithered data may prefer coarser values. +""" + +from .adapter import InstrumentAdapter, register + +WFC3_IR = register( + InstrumentAdapter( + key="wfc3_ir", + mast_instrument_name="WFC3/IR", + native_scale=0.128, + calibrated_suffix="FLT", + reference_env_key="iref", + crds_reference_subpath="references/hst/wfc3", + supports_cte_correction=False, + default_drizzle_kwargs={ + "skymethod": "globalmin+match", + "final_wht_type": "IVM", + "final_units": "cps", + "final_rot": 0.0, + }, + saturation_dn=78_000.0, + recommended_final_scale=0.065, + ) +) diff --git a/autoreduce/instruments/wfc3_uvis.py b/autoreduce/instruments/wfc3_uvis.py new file mode 100644 index 0000000..2596507 --- /dev/null +++ b/autoreduce/instruments/wfc3_uvis.py @@ -0,0 +1,30 @@ +""" +WFC3/UVIS — adapter #2 (roadmap phase 2). + +The ACS-like path: CTE-corrected ``_flc`` exposures, ``iref`` references, +cps/IVM/north-up outputs. Native plate scale 0.0396″/pix — the published +lensing anchor is the Bayer et al. (arXiv:1803.05952) F390W reduction of +SDSS J0252+0039 at exactly this output scale with pixfrac 1.0. +""" + +from .adapter import InstrumentAdapter, register + +WFC3_UVIS = register( + InstrumentAdapter( + key="wfc3_uvis", + mast_instrument_name="WFC3/UVIS", + native_scale=0.0396, + calibrated_suffix="FLC", + reference_env_key="iref", + crds_reference_subpath="references/hst/wfc3", + supports_cte_correction=True, + default_drizzle_kwargs={ + "skymethod": "globalmin+match", + "final_wht_type": "IVM", + "final_units": "cps", + "final_rot": 0.0, + }, + saturation_dn=63_000.0, + recommended_final_scale=0.0396, + ) +) diff --git a/autoreduce/pipeline.py b/autoreduce/pipeline.py index 00180f0..ad2f6fc 100644 --- a/autoreduce/pipeline.py +++ b/autoreduce/pipeline.py @@ -118,11 +118,21 @@ def reduce_target( target_xy = WCS(header).world_to_pixel_values(spec.ra, spec.dec) selection = stars_mod.StarSelection() + # Saturation is per exposure, not per stack: a star saturates when its + # rate fills the well within one exposure, so the cps cap divides the + # full well by the longest single-exposure time — never the mosaic total. + max_single_exptime = max( + float(fits.getheader(p).get("EXPTIME", 0.0)) for p in exposures + ) + if max_single_exptime <= 0.0: + raise ValueError("no exposure carries a positive EXPTIME header") stars = stars_mod.find_stars( sci, selection, target_xy=(float(target_xy[0]), float(target_xy[1])), - peak_max=selection.saturation_fraction * adapter.saturation_dn / float(exptime), + peak_max=selection.saturation_fraction + * adapter.saturation_dn + / max_single_exptime, ) psf, psf_full, psf_diag = epsf_mod.build_epsf( sci, stars, spec.psf_shape, spec.psf_full_shape diff --git a/docs/design/roadmap.md b/docs/design/roadmap.md index 06769ca..deb9a91 100644 --- a/docs/design/roadmap.md +++ b/docs/design/roadmap.md @@ -13,17 +13,14 @@ geometry, calibrated-product naming (`_flc` vs `_flt` vs `_cal`), units and gain, recommended combine parameters, and the PSF-model source. ACS/WFC is adapter #1; nothing outside `instruments/` may mention a detector by name. -## HST/WFC3 (IR + UVIS) - -- Same stage graph; new adapter. UVIS is ACS-like (CTE-corrected `_flc`, - drizzle); IR differs (no CTE correction, up-the-ramp sampling, different - DQ semantics, 0.13″ native pixels → different `final_scale` choice). -- The [ajshajib/hst-lens](https://github.com/ajshajib/hst-lens) notebooks - (download / IR / UVIS) serve as a step checklist to audit against — not as - architecture. Their gaps (no ACS, notebook-form, unclear noise treatment) - are what this project exists to fix. -- Other ACS filters (F435W, F606W…) are config, not code: the adapter already - parameterizes filter-dependent zero-points and PSF models. +## HST/WFC3 (IR + UVIS) — **in progress (phase 2, PyAutoReduce#4)** + +- Design deltas live in [`wfc3.md`](wfc3.md); adapters `wfc3_uvis` / + `wfc3_ir` implemented. UVIS is ACS-like (CTE-corrected `_flc`); IR differs + (`_flt`, no CTE correction, up-the-ramp CR rejection, 0.128″ native → + recommended 0.065″ output). +- Other ACS/WFC3 filters (F435W, F606W…) are config, not code: the adapter + already parameterizes the filter-dependent pieces. ## JWST (NIRCam first) diff --git a/docs/design/wfc3.md b/docs/design/wfc3.md new file mode 100644 index 0000000..3a658b7 --- /dev/null +++ b/docs/design/wfc3.md @@ -0,0 +1,65 @@ +# WFC3 (UVIS + IR) — per-stage deltas vs the ACS/WFC design + +Phase 2. `hst_acs_pipeline.md` remains the stage-by-stage spec; this page +records only what WFC3 changes. Everything below lives in the two adapters +(`instruments/wfc3_uvis.py`, `instruments/wfc3_ir.py`) — no stage module +mentions a detector. + +## WFC3/UVIS — the ACS-like path + +| Stage | Delta vs ACS | +|-------|--------------| +| acquire | `_flc` (CTE-corrected) as ACS; references key **`iref`**, CRDS subpath `references/hst/wfc3` | +| align / drizzle / noise / psf / package | unchanged — same recipes, same dials | +| scale | native **0.0396″/pix**; adapter recommends output at native scale | +| saturation | ~63 ke- full well (star selection peak cut) | + +**Validation anchor:** the published [Bayer et al. +(arXiv:1803.05952)](https://arxiv.org/abs/1803.05952) F390W reduction of +SDSS J0252+0039 — output 0.0396″/pix, **pixfrac 1.0**, noise recipe +σ = √(N/W + σ²_sky) with σ_sky ≈ 0.002 e-/s. Our integration script +(`scripts/reduce_j0252_wfc3.py --channel uvis`) reduces the same data with those dials and +checks: output units e-/s, empirical σ_sky in that regime, noise-map +consistency, WHT uniformity. (Note their correlated-noise treatment is +blank-sky *realizations*, not the scalar R — with pixfrac 1.0 at native +scale, R = 1.5; comparisons account for whether R is applied.) + +## WFC3/IR — the genuinely different path + +| Stage | Delta vs ACS | +|-------|--------------| +| acquire | **`_flt`** — no CTE correction exists for the IR channel; `iref` references | +| drizzle | up-the-ramp fitting in `calwf3` already rejects most CRs per read; `driz_cr` still runs on multi-exposure stacks for the residue (defaults-first) — documented, revisit if IR integrations show over-flagging | +| scale | native **0.128″/pix** under-samples the PSF; adapter recommends **0.065″/pix** for dithered data (half-native, in the 0.06–0.08 deep-field range). The dial stays user-facing; the fine-grid Casertano branch (s < p) then applies, so R is materially larger — reported per run as always | +| units | `_flt` IR data are already e-/s (count rates); `final_units='cps'` unchanged | +| saturation | ~78 ke- effective full well | +| psf | same tiers; STScI focus-diverse ePSF grids exist for IR when tier 2 lands | + +## Coverage audit vs `ajshajib/hst-lens` (the checklist, not the architecture) + +Their three notebooks (Download / IR / UVIS) cover: archive download, +per-channel calibrated products, AstroDrizzle combination, and cutouts. Ours +adds what they lack: instrument adapters (theirs is notebook-per-channel), +provenance (`reduction.json`), an explicit noise-map recipe with correlated- +noise handling, tiered PSF construction with diagnostics, cache/eviction, and +loud-failure contracts. Nothing in their steps is absent from our stage graph. + +## Integration finding — IR pixfrac/coverage (2026-07-08) + +Reducing the J0252+0039 F160W snapshot (program 11202) at the recommended +0.065″/pix with the phase-1 default pixfrac 0.8 left **230 zero-weight +speckle pixels inside the cutout** — the few-dither + fine-grid (s ≈ 0.51 < +p) regime — and the finite-noise packaging guard refused to ship the dataset, +exactly as designed. pixfrac 1.0 closes coverage at the cost of a larger +correlated-noise factor (reported per run, as always). Rule of thumb recorded +here: **on the IR channel, few-dither data at sub-native output scales needs +pixfrac → 1.0 (or a coarser scale dial)**; the WHT-uniformity diagnostic and +the finite-noise guard enforce the trade-off loudly rather than letting a +holey dataset through. + +## Open items + +- IR integration target: discovered via MAST at run time; if no suitable + IR lens dataset is reachable, the leg parks as a batched question. +- Tier-2 PSF for IR (focus-diverse ePSF grids) — with the roadmap's tier-2 + work, not phase 2. diff --git a/scripts/reduce_j0252_wfc3.py b/scripts/reduce_j0252_wfc3.py new file mode 100644 index 0000000..0ac5799 --- /dev/null +++ b/scripts/reduce_j0252_wfc3.py @@ -0,0 +1,147 @@ +""" +WFC3 integration + acceptance (issue #4): SDSS J0252+0039. + +--channel uvis : F390W through wfc3_uvis at the published Bayer et al. + (arXiv:1803.05952) dials — 0.0396"/pix, pixfrac 1.0 — then + check our products against their published numbers (units + e-/s; sigma_sky ~ 0.002 e-/s; noise closure; R accounting). +--channel ir : discover WFC3/IR imaging at the same coordinates via MAST + and reduce it at the adapter-recommended 0.065"/pix; + internal validation only (no published anchor). + +Run: ~/venv/PyAuto/bin/python scripts/reduce_j0252_wfc3.py --channel uvis +Network + drizzlepac required; unit tests never import this. +""" + +import argparse +import json +import sys +from pathlib import Path + +import numpy as np + +REPO = Path(__file__).resolve().parent.parent +sys.path.insert(0, str(REPO)) + +from autoreduce import TargetSpec, reduce_target # noqa: E402 +from autoreduce import instruments # noqa: E402 + +# SDSS J0252+0039: 02h52m45.21s +00d39m58.4s (legacy dataset info.json) +RA, DEC = 43.188375, 0.666222 +CACHE_ROOT = REPO / "scripts" / "cache" +OUTPUT_ROOT = REPO / "scripts" / "output" + +BAYER_SIGMA_SKY = 0.002 # e-/s, published for the F390W reduction + + +def spec_for(channel: str) -> TargetSpec: + if channel == "uvis": + return TargetSpec( + name="j0252+0039_f390w", + ra=RA, + dec=DEC, + instrument="wfc3_uvis", + filter_name="F390W", + final_scale=0.0396, # Bayer dial + final_pixfrac=1.0, # Bayer dial + ) + if channel == "ir": + filter_name = discover_ir_filter() + return TargetSpec( + name=f"j0252+0039_{filter_name.lower()}", + ra=RA, + dec=DEC, + instrument="wfc3_ir", + filter_name=filter_name, + final_scale=instruments.get("wfc3_ir").recommended_final_scale, + # Few-dither snapshot on a half-native grid: pixfrac 0.8 leaves + # zero-weight speckle (the finite-noise guard caught it); the + # full drop closes coverage at the cost of a larger R — the + # dial trade-off working as documented. + final_pixfrac=1.0, + # 281 px at 0.065" spans 18.3" — far more sky than the ACS-era + # 14" footprint, and it clips a zero-coverage detector-defect + # blob 8.5" from the lens (the guard refused it). Match the ACS + # sky footprint instead: 14" / 0.065 -> 215 px. + cutout_shape=(215, 215), + ) + raise ValueError(channel) + + +def discover_ir_filter() -> str: + """Find which WFC3/IR filter (if any) covers the target.""" + from astropy.coordinates import SkyCoord + from astroquery.mast import Observations + + from autoreduce.acquire.mast import select_observations + + obs = Observations.query_criteria( + coordinates=SkyCoord(RA, DEC, unit="deg"), + radius="0.5 arcmin", + obs_collection="HST", + instrument_name="WFC3/IR", + dataproduct_type="image", + ) + direct = select_observations(obs) + if not direct: + sys.exit( + "[ir] no direct WFC3/IR observations at J0252+0039 — the IR leg " + "needs a different target (parked question on issue #4)" + ) + # HAP composite rows carry the pseudo-filter 'detection'; never reduce it. + filters = sorted( + {str(row["filters"]) for row in direct} - {"detection"} + ) + if not filters: + sys.exit("[ir] only HAP composite products found — no real filter") + print(f"[ir] direct WFC3/IR observations found; filters: {filters}") + preferred = [f for f in filters if f == "F160W"] or filters + return preferred[0] + + +def validate(channel: str, record: dict, out_dir: Path): + from astropy.io import fits + + noise = fits.getdata(out_dir / "noise_map.fits").astype(float) + data = fits.getdata(out_dir / "data.fits").astype(float) + with fits.open(out_dir / "data.fits") as hdul: + bunit = hdul[0].header.get("BUNIT", "unknown") + + r_factor = record["noise"]["correlated_noise_factor"] + sky_rms = record["noise"]["empirical_background_rms"] + summary = { + "channel": channel, + "n_exposures": record["acquire"]["n_exposures"], + "bunit": bunit, + "weight_uniformity": record["drizzle"]["weight_uniformity"], + "correlated_noise_factor": r_factor, + "empirical_sky_rms_cps": sky_rms, + "noise_map_min_median": [float(np.nanmin(noise)), float(np.nanmedian(noise))], + "psf": record["psf"], + } + if channel == "uvis": + # Published anchor: Bayer et al. sigma_sky ~ 0.002 e-/s at these dials. + # Our empirical sky RMS is pre-R; theirs enters sigma_n pre-realization. + summary["bayer_sigma_sky_cps"] = BAYER_SIGMA_SKY + summary["sky_rms_over_bayer"] = sky_rms / BAYER_SIGMA_SKY + # Noise-floor closure: in blank sky the map should approach R * sky_rms. + summary["noise_floor_over_R_times_sky"] = float( + np.nanpercentile(noise, 5) / (r_factor * sky_rms) + ) + print(f"[{channel}] ---- validation ----") + print(json.dumps(summary, indent=2)) + (out_dir / "validation_summary.json").write_text(json.dumps(summary, indent=2)) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--channel", required=True, choices=["uvis", "ir"]) + args = parser.parse_args() + + spec = spec_for(args.channel) + record = reduce_target(spec, cache_root=CACHE_ROOT, output_root=OUTPUT_ROOT) + validate(args.channel, record, OUTPUT_ROOT / spec.name) + + +if __name__ == "__main__": + main() diff --git a/test_autoreduce/test_wfc3.py b/test_autoreduce/test_wfc3.py new file mode 100644 index 0000000..a68d100 --- /dev/null +++ b/test_autoreduce/test_wfc3.py @@ -0,0 +1,79 @@ +import pytest + +from autoreduce import instruments +from autoreduce.drizzle.combine import drizzle_kwargs_for +from autoreduce.target import TargetSpec + + +class TestWFC3Adapters: + def test_both_channels_registered(self): + assert "wfc3_uvis" in instruments.registered_keys() + assert "wfc3_ir" in instruments.registered_keys() + + def test_uvis_is_acs_like(self): + uvis = instruments.get("wfc3_uvis") + assert uvis.calibrated_suffix == "FLC" + assert uvis.supports_cte_correction + assert uvis.reference_env_key == "iref" + assert uvis.native_scale == pytest.approx(0.0396) + assert uvis.recommended_final_scale == pytest.approx(0.0396) + + def test_ir_is_the_different_path(self): + ir = instruments.get("wfc3_ir") + # No CTE correction exists for the IR channel: _flt, not _flc. + assert ir.calibrated_suffix == "FLT" + assert not ir.supports_cte_correction + assert ir.reference_env_key == "iref" + assert ir.native_scale == pytest.approx(0.128) + # Under-sampled detector: recommendation is a finer output grid. + assert ir.recommended_final_scale < ir.native_scale + + def test_wfc3_crds_subpath_owned_by_adapter(self): + for key in ("wfc3_uvis", "wfc3_ir"): + assert instruments.get(key).crds_reference_subpath == "references/hst/wfc3" + + def test_acs_recommendation_unchanged(self): + # Phase-1 regression: the ACS path must not change. + assert instruments.get("acs_wfc").recommended_final_scale == pytest.approx(0.05) + + +class TestWFC3DrizzleKwargs: + def test_uvis_kwargs_at_bayer_scale(self): + spec = TargetSpec( + name="j0252", ra=43.19, dec=0.666, instrument="wfc3_uvis", + filter_name="F390W", final_scale=0.0396, final_pixfrac=1.0, + ) + kwargs = drizzle_kwargs_for(spec, instruments.get("wfc3_uvis"), 4) + assert kwargs["final_scale"] == pytest.approx(0.0396) + assert kwargs["final_pixfrac"] == pytest.approx(1.0) + assert kwargs["final_units"] == "cps" + assert kwargs["driz_cr"] + + def test_ir_single_exposure_branch_still_applies(self): + spec = TargetSpec( + name="x", ra=0.0, dec=0.0, instrument="wfc3_ir", + filter_name="F160W", final_scale=0.065, + ) + kwargs = drizzle_kwargs_for(spec, instruments.get("wfc3_ir"), 1) + assert not kwargs["driz_cr"] + + def test_ir_fine_grid_correlation_factor(self): + # Drizzling 0.128 -> 0.065 with pixfrac 0.8: s < p branch engaged. + from autoreduce.noise.rms import casertano_r + + ir = instruments.get("wfc3_ir") + s = ir.scale_ratio(0.065) + assert s < 0.8 + assert casertano_r(0.8, s) > casertano_r(0.8, 1.0) + + +def test_crds_environment_uses_adapter_subpath(tmp_path, monkeypatch): + import os + + from autoreduce.acquire.crds import configure_environment + + monkeypatch.delenv("iref", raising=False) + monkeypatch.delenv("CRDS_PATH", raising=False) + env = configure_environment(tmp_path, instruments.get("wfc3_ir")) + assert env["iref"].endswith("references/hst/wfc3/") + assert os.environ["iref"] == env["iref"]