diff --git a/autoreduce/acquire/crds.py b/autoreduce/acquire/crds.py index 1db08e6..0e0b2d5 100644 --- a/autoreduce/acquire/crds.py +++ b/autoreduce/acquire/crds.py @@ -15,26 +15,28 @@ from ..instruments import InstrumentAdapter -CRDS_SERVER_URL = "https://hst-crds.stsci.edu" - def configure_environment(references_root: Path, adapter: InstrumentAdapter) -> dict: """ - Set the CRDS variables for this process. Must run before drizzlepac is - imported anywhere in the process. Returns the mapping applied. + Set the CRDS variables for this process. Must run before drizzlepac / + the jwst pipeline are imported anywhere in the process. Returns the + mapping applied. Deliberately overrides any inherited CRDS_PATH/jref: the pipeline is a pure function of the target spec plus the archive, so its reference files live in *its* cache, not wherever the shell environment happens to point. + The server URL is the adapter's (hst-crds vs jwst-crds). """ env = { - "CRDS_SERVER_URL": CRDS_SERVER_URL, + "CRDS_SERVER_URL": adapter.crds_server_url, "CRDS_PATH": str(references_root), - adapter.reference_env_key: str( - Path(references_root) / adapter.crds_reference_subpath - ) - + "/", } + # HST-style tools resolve references through an iraf-style variable + # (jref$/iref$); the jwst pipeline reads CRDS_PATH directly. + if adapter.reference_env_key != "CRDS_PATH": + env[adapter.reference_env_key] = ( + str(Path(references_root) / adapter.crds_reference_subpath) + "/" + ) os.environ.update(env) return env diff --git a/autoreduce/acquire/footprint.py b/autoreduce/acquire/footprint.py new file mode 100644 index 0000000..32cbcad --- /dev/null +++ b/autoreduce/acquire/footprint.py @@ -0,0 +1,63 @@ +""" +Detector-footprint filtering (design docs stage 1; added for JWST). + +A survey visit's exposures span many detectors, most of which never touch the +target; combining them wastes memory and time (and on this machine, OOMs the +jwst pipeline). Keep only calibrated exposures whose detector footprint +contains the target, with a margin for the cutout and dither pattern. + +Uses the approximate FITS WCS every calibrated product carries in its SCI +extension (JWST cal files carry both gwcs and FITS-approx; HST flt/flc carry +FITS WCS) — footprint containment at arcsecond precision, which is all this +filter needs. +""" + +from pathlib import Path +from typing import List, Tuple + + +def covers_target(path: Path, ra: float, dec: float, margin_arcsec: float) -> bool: + """True if any SCI extension's footprint contains (ra, dec) ± margin. + + Containment is tested in *pixel* space (project the target through the + extension's WCS, allow a margin in pixels) — immune to the RA-wraparound + and cos(dec) pitfalls of sky-coordinate bounding boxes. + """ + import numpy as np + from astropy.io import fits + from astropy.wcs import WCS + from astropy.wcs.utils import proj_plane_pixel_scales + + with fits.open(path) as hdul: + for hdu in hdul: + if hdu.name != "SCI" or hdu.data is None: + continue + wcs = WCS(hdu.header, naxis=2) + ny, nx = hdu.data.shape[-2], hdu.data.shape[-1] + x, y = wcs.world_to_pixel_values(ra, dec) + if not (np.isfinite(x) and np.isfinite(y)): + continue + scale_arcsec = float( + np.mean(proj_plane_pixel_scales(wcs)) * 3600.0 + ) + m = margin_arcsec / scale_arcsec + if -m <= float(x) <= nx - 1 + m and -m <= float(y) <= ny - 1 + m: + return True + return False + + +def filter_to_target( + exposures: List[Path], ra: float, dec: float, margin_arcsec: float = 30.0 +) -> Tuple[List[Path], List[Path]]: + """Split exposures into (covering, skipped); loud if nothing covers.""" + covering, skipped = [], [] + for path in exposures: + (covering if covers_target(path, ra, dec, margin_arcsec) else skipped).append( + path + ) + if not covering: + raise LookupError( + f"none of {len(exposures)} exposures cover ({ra}, {dec}) — " + f"wrong coordinates or a query/footprint bug" + ) + return covering, skipped diff --git a/autoreduce/acquire/mast.py b/autoreduce/acquire/mast.py index 1fc694a..8f56499 100644 --- a/autoreduce/acquire/mast.py +++ b/autoreduce/acquire/mast.py @@ -56,7 +56,7 @@ def query_exposures( obs = Observations.query_criteria( coordinates=coord, radius=radius, - obs_collection="HST", + obs_collection=adapter.mast_obs_collection, instrument_name=adapter.mast_instrument_name, filters=filter_name, dataproduct_type="image", diff --git a/autoreduce/drizzle/combine.py b/autoreduce/drizzle/combine.py index ce687eb..fef6b42 100644 --- a/autoreduce/drizzle/combine.py +++ b/autoreduce/drizzle/combine.py @@ -48,10 +48,21 @@ def combine( output_dir: Path, ) -> Tuple[Path, Path, Dict]: """ - Run AstroDrizzle; return (sci_path, wht_path, provenance_fragment). + Combine exposures via the adapter's backend; return + (sci_path, wht_path, provenance_fragment). Requires the CRDS environment configured (acquire.crds) beforehand. """ + if adapter.combine_backend == "jwst_image3": + from . import jwst_combine + + return jwst_combine.combine(exposures, spec, adapter, output_dir) + if adapter.combine_backend != "astrodrizzle": + raise ValueError( + f"unknown combine backend {adapter.combine_backend!r} " + f"for instrument {adapter.key}" + ) + from astropy.io import fits from drizzlepac import astrodrizzle diff --git a/autoreduce/drizzle/jwst_combine.py b/autoreduce/drizzle/jwst_combine.py new file mode 100644 index 0000000..589d047 --- /dev/null +++ b/autoreduce/drizzle/jwst_combine.py @@ -0,0 +1,113 @@ +""" +JWST combination backend (roadmap phase 3): calwebb_image3 — the drizzle +analogue (tweakreg / skymatch / outlier_detection / resample). Defaults-first: +the pipeline runs with its own defaults; only the lensing dials map through +(output pixel scale, pixfrac, kernel, north-up rotation, IVM weighting). + +The ``_i2d`` product is multi-extension (SCI/ERR/CON/WHT/VAR_*); this module +normalizes it to the package's internal contract — standalone ``sci``/``wht``/ +``err`` FITS files with the WCS and an EXPTIME key — so every downstream +stage (noise, psf, package) is backend-agnostic. +""" + +import json +from pathlib import Path +from typing import Dict, List, Tuple + +import numpy as np + +from ..instruments import InstrumentAdapter +from ..target import TargetSpec +from .diagnostics import check_weight_uniformity +from ..noise.rms import casertano_r + + +def combine( + exposures: List[Path], + spec: TargetSpec, + adapter: InstrumentAdapter, + output_dir: Path, +) -> Tuple[Path, Path, Dict]: + """Run calwebb_image3; return (sci_path, wht_path, provenance_fragment).""" + from astropy.io import fits + from jwst.associations.asn_from_list import asn_from_list + from jwst.associations.lib.rules_level3_base import DMS_Level3_Base + from jwst.pipeline import Image3Pipeline + + # Resolve before the chdir below: relative paths would dangle afterwards. + output_dir = Path(output_dir).resolve() + output_dir.mkdir(parents=True, exist_ok=True) + product_name = f"{spec.name}_{spec.filter_name}".lower() + + asn = asn_from_list( + [str(p) for p in exposures], rule=DMS_Level3_Base, product_name=product_name + ) + asn_path = output_dir / f"{product_name}_asn.json" + _, serialized = asn.dump(format="json") + asn_path.write_text(serialized) + + import os + + cwd = os.getcwd() + os.chdir(output_dir) + try: + Image3Pipeline.call( + str(asn_path), + output_dir=str(output_dir), + save_results=True, # stpipe .call() discards results otherwise + in_memory=False, # on-disk models: image3 OOMs this machine otherwise + steps={ + "resample": { + "pixel_scale": spec.final_scale, + "pixfrac": spec.final_pixfrac, + "kernel": spec.final_kernel, + "rotation": 0.0, + "weight_type": adapter.default_drizzle_kwargs.get( + "weight_type", "ivm" + ), + }, + }, + ) + finally: + os.chdir(cwd) + + i2d = output_dir / f"{product_name}_i2d.fits" + if not i2d.exists(): + raise FileNotFoundError(f"calwebb_image3 did not produce {i2d}") + + # Normalize to the internal contract: standalone sci/wht/err files. + with fits.open(i2d) as hdul: + sci = hdul["SCI"].data.astype(np.float32) + err = hdul["ERR"].data.astype(np.float32) + wht = hdul["WHT"].data.astype(np.float32) + header = hdul["SCI"].header.copy() + header["EXPTIME"] = hdul[0].header.get( + "XPOSURE", hdul["SCI"].header.get("XPOSURE", 0.0) + ) + header["BUNIT"] = hdul["SCI"].header.get("BUNIT", "MJy/sr") + + paths = {} + for name, data in (("sci", sci), ("wht", wht), ("err", err)): + path = output_dir / f"{product_name}_{name}.fits" + fits.PrimaryHDU(data, header=header).writeto(path, overwrite=True) + paths[name] = path + + provenance = { + "backend": "jwst_image3", + "n_exposures": len(exposures), + "exposures": [Path(p).name for p in exposures], + "single_exposure_branch": len(exposures) == 1, + "resample_kwargs": { + "pixel_scale": spec.final_scale, + "pixfrac": spec.final_pixfrac, + "kernel": spec.final_kernel, + "rotation": 0.0, + }, + "correlated_noise_factor": casertano_r( + spec.final_pixfrac, adapter.scale_ratio(spec.final_scale) + ), + "weight_uniformity": check_weight_uniformity(wht), + "err_path": str(paths["err"]), + "i2d_path": str(i2d), + } + return paths["sci"], paths["wht"], provenance diff --git a/autoreduce/instruments/__init__.py b/autoreduce/instruments/__init__.py index 2cd02f7..7ea8b65 100644 --- a/autoreduce/instruments/__init__.py +++ b/autoreduce/instruments/__init__.py @@ -9,3 +9,4 @@ from .acs_wfc import ACS_WFC from .wfc3_uvis import WFC3_UVIS from .wfc3_ir import WFC3_IR +from .nircam import NIRCAM_SW, NIRCAM_LW, nircam_adapter_for_filter diff --git a/autoreduce/instruments/adapter.py b/autoreduce/instruments/adapter.py index 8bcd37b..5a63c62 100644 --- a/autoreduce/instruments/adapter.py +++ b/autoreduce/instruments/adapter.py @@ -24,6 +24,12 @@ class InstrumentAdapter: # 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 + # Observatory-level routing (phase 3): which archive/CRDS ecosystem and + # which combination backend this instrument reduces through. + observatory: str = "hst" # "hst" | "jwst" + crds_server_url: str = "https://hst-crds.stsci.edu" + combine_backend: str = "astrodrizzle" # "astrodrizzle" | "jwst_image3" + mast_obs_collection: str = "HST" def scale_ratio(self, final_scale: float) -> float: """s = output scale / native scale, as used by the Casertano factor.""" diff --git a/autoreduce/instruments/nircam.py b/autoreduce/instruments/nircam.py new file mode 100644 index 0000000..8fb9b4f --- /dev/null +++ b/autoreduce/instruments/nircam.py @@ -0,0 +1,70 @@ +""" +JWST/NIRCam — adapters #4 and #5 (roadmap phase 3), the first non-HST path. + +Both channels reduce from level-2 ``_cal`` products (calwebb_image2 output, +MJy/sr) and combine through the ``jwst`` pipeline's calwebb_image3 — the +drizzle analogue (tweakreg / skymatch / outlier_detection / resample). CRDS +routes through the JWST server; references sync under ``references/jwst``. + +Output-scale recommendations follow the COSMOS-Web mosaic convention the +parity anchor uses: SW at 0.03″/pix, LW at 0.06″/pix. ``saturation_dn`` is +the detector full well in electrons, as elsewhere; NIRCam ~ 105 ke- (2RG +arrays), conservative here. +""" + +from .adapter import InstrumentAdapter, register + +_COMMON = dict( + calibrated_suffix="CAL", + reference_env_key="CRDS_PATH", # jwst pipeline reads CRDS_PATH directly + crds_reference_subpath="references/jwst", + supports_cte_correction=False, + observatory="jwst", + crds_server_url="https://jwst-crds.stsci.edu", + combine_backend="jwst_image3", + mast_obs_collection="JWST", + default_drizzle_kwargs={ + # calwebb_image3 resample step keywords (dial-mapped in jwst_combine) + "pixfrac": 1.0, + "kernel": "square", + "weight_type": "ivm", + }, + saturation_dn=100_000.0, +) + +NIRCAM_SW = register( + InstrumentAdapter( + key="nircam_sw", + mast_instrument_name="NIRCAM/IMAGE", + native_scale=0.031, + recommended_final_scale=0.03, + **_COMMON, + ) +) + +NIRCAM_LW = register( + InstrumentAdapter( + key="nircam_lw", + mast_instrument_name="NIRCAM/IMAGE", + native_scale=0.063, + recommended_final_scale=0.06, + **_COMMON, + ) +) + +# NIRCam filter -> channel routing (wavelength < 2.4 micron = SW). +SW_FILTERS = {"F070W", "F090W", "F115W", "F140M", "F150W", "F162M", "F164N", + "F150W2", "F182M", "F187N", "F200W", "F210M", "F212N"} +LW_FILTERS = {"F250M", "F277W", "F300M", "F322W2", "F323N", "F335M", "F356W", + "F360M", "F405N", "F410M", "F430M", "F444W", "F460M", "F466N", + "F470N", "F480M"} + + +def nircam_adapter_for_filter(filter_name: str) -> InstrumentAdapter: + """Route a NIRCam filter to its channel adapter; loud on unknown filters.""" + name = filter_name.upper() + if name in SW_FILTERS: + return NIRCAM_SW + if name in LW_FILTERS: + return NIRCAM_LW + raise KeyError(f"unknown NIRCam filter {filter_name!r} — not in SW/LW tables") diff --git a/autoreduce/noise/jwst_rms.py b/autoreduce/noise/jwst_rms.py new file mode 100644 index 0000000..17fe0fe --- /dev/null +++ b/autoreduce/noise/jwst_rms.py @@ -0,0 +1,48 @@ +""" +JWST noise stage (roadmap phase 3): *read, don't construct*. + +The resampled ``_i2d`` product already carries a per-pixel total-error map +(ERR: Poisson + read noise + flat, propagated and resampled by the jwst +pipeline), so stage 4 reads it rather than rebuilding it from weights — and +then applies the same correlated-noise factor R the HST path uses, since +resample correlates neighbouring pixels exactly as drizzle does and the +propagated ERR is a per-pixel quantity that underestimates the effective +noise a lens-model chi^2 sees. + +A consistency check against the empirical blank-sky RMS of the mosaic is +recorded in provenance; large disagreement means the upstream error model +and the sky disagree and must be investigated, not absorbed. +""" + +from pathlib import Path +from typing import Dict, Tuple + +import numpy as np + +from .rms import empirical_background_rms + + +def noise_map_from_error( + err: np.ndarray, + sci: np.ndarray, + correlated_noise_factor: float = 1.0, +) -> Tuple[np.ndarray, Dict]: + """RMS map from a propagated ERR array; NaN/zero stay NaN (loud later).""" + if err.shape != sci.shape: + raise ValueError(f"shape mismatch: err {err.shape} vs sci {sci.shape}") + if correlated_noise_factor < 1.0: + raise ValueError( + f"correlated-noise factor must be >= 1: {correlated_noise_factor}" + ) + noise = np.where( + np.isfinite(err) & (err > 0.0), correlated_noise_factor * err, np.nan + ) + + sky_rms = empirical_background_rms(sci[np.isfinite(noise)]) + err_floor = float(np.nanpercentile(noise, 5)) / correlated_noise_factor + consistency = { + "empirical_sky_rms": sky_rms, + "err_5th_percentile_pre_R": err_floor, + "sky_over_err_floor": sky_rms / err_floor if err_floor > 0 else float("inf"), + } + return noise, consistency diff --git a/autoreduce/noise/rms.py b/autoreduce/noise/rms.py index 25446a2..a213d6f 100644 --- a/autoreduce/noise/rms.py +++ b/autoreduce/noise/rms.py @@ -77,6 +77,75 @@ def assert_finite_within(noise_map: np.ndarray, region_name: str) -> None: ) +# Masked-by-noise convention: bad pixels carry effectively infinite noise so +# any chi^2 ignores them — the same treatment the legacy noise-scaled +# datasets use for artifacts and contaminants. +MASKED_NOISE_VALUE = 1.0e8 + + +def mask_isolated_bad_pixels( + data_cut: np.ndarray, + noise_cut: np.ndarray, + center_xy, + pixel_scale: float, + max_bad_fraction: float = 0.005, + protect_radius_arcsec: float = 1.5, + region_name: str = "cutout", +): + """ + Apply the bad-pixel policy to a cutout pair. + + Isolated non-finite/non-positive noise pixels (fully-rejected or dead + pixels — routine in deep resampled stacks) are set to `MASKED_NOISE_VALUE` + with the data zeroed, and the count/positions are returned for provenance. + The failure stays loud where it matters: more than `max_bad_fraction` of + the cutout bad, or any bad pixel within `protect_radius_arcsec` of the + target centre (the lens itself must be clean). + """ + bad = ~np.isfinite(noise_cut) | (noise_cut <= 0.0) + n_bad = int(bad.sum()) + if n_bad == 0: + return data_cut, noise_cut, {"n_masked_pixels": 0} + + # "Isolated" is enforced: a bad pixel with two or more bad 4-neighbours + # marks a structured defect (blob/column), which must fail loudly — only + # scattered singletons and pairs are maskable. + neighbours = sum( + np.roll(bad, shift, axis) for shift, axis in ((1, 0), (-1, 0), (1, 1), (-1, 1)) + ) + if (bad & (neighbours >= 2)).any(): + raise ValueError( + f"structured bad-pixel region in {region_name} ({n_bad} bad px with " + f"contiguous clustering) — fix the reduction, don't mask a defect" + ) + + fraction = n_bad / bad.size + if fraction > max_bad_fraction: + raise ValueError( + f"{n_bad} bad noise pixels ({fraction:.2%}) in {region_name} exceed " + f"the {max_bad_fraction:.2%} policy limit — fix the reduction " + f"(coverage, weights), don't mask wholesale" + ) + ys, xs = np.where(bad) + cx, cy = center_xy + r_arcsec = np.hypot(ys - cy, xs - cx) * pixel_scale + if (r_arcsec < protect_radius_arcsec).any(): + raise ValueError( + f"bad noise pixel within {protect_radius_arcsec}\" of the target " + f"centre in {region_name} — the lens region must reduce cleanly" + ) + + data_out = np.where(bad, 0.0, data_cut) + noise_out = np.where(bad, MASKED_NOISE_VALUE, noise_cut) + diagnostics = { + "n_masked_pixels": n_bad, + "masked_fraction": fraction, + "masked_noise_value": MASKED_NOISE_VALUE, + "min_masked_radius_arcsec": float(r_arcsec.min()), + } + return data_out, noise_out, diagnostics + + 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/package/cutout.py b/autoreduce/package/cutout.py index 4436446..f02cd38 100644 --- a/autoreduce/package/cutout.py +++ b/autoreduce/package/cutout.py @@ -11,18 +11,15 @@ import numpy as np -def cutout_to_fits( +def make_cutout( data: np.ndarray, header, ra: float, dec: float, shape: Tuple[int, int], - out_path: Path, - extra_header: dict = None, -) -> np.ndarray: - """Cut `shape` around (ra, dec) and write with an intact cutout WCS.""" +): + """Cut `shape` around (ra, dec); returns (cut_data, out_header, center_xy).""" from astropy.coordinates import SkyCoord - from astropy.io import fits from astropy.nddata import Cutout2D from astropy.wcs import WCS @@ -33,11 +30,31 @@ def cutout_to_fits( for key in ("BUNIT", "EXPTIME", "TEXPTIME", "FILTER", "INSTRUME", "TELESCOP"): if key in header: out_header[key] = header[key] - if extra_header: - for key, value in extra_header.items(): - out_header[key] = value + center = cut.wcs.world_to_pixel_values(ra, dec) + return cut.data, out_header, (float(center[0]), float(center[1])) + + +def write_fits(data: np.ndarray, header, out_path: Path) -> None: + from astropy.io import fits - fits.PrimaryHDU(cut.data.astype(np.float32), header=out_header).writeto( + fits.PrimaryHDU(data.astype(np.float32), header=header).writeto( out_path, overwrite=True ) - return cut.data + + +def cutout_to_fits( + data: np.ndarray, + header, + ra: float, + dec: float, + shape: Tuple[int, int], + out_path: Path, + extra_header: dict = None, +) -> np.ndarray: + """Cut and write in one step (kept for callers with no post-processing).""" + cut_data, out_header, _ = make_cutout(data, header, ra, dec, shape) + if extra_header: + for key, value in extra_header.items(): + out_header[key] = value + write_fits(cut_data, out_header, out_path) + return cut_data diff --git a/autoreduce/pipeline.py b/autoreduce/pipeline.py index ad2f6fc..444d6e1 100644 --- a/autoreduce/pipeline.py +++ b/autoreduce/pipeline.py @@ -16,6 +16,7 @@ from . import instruments from .acquire import cache as cache_mod from .acquire import crds as crds_mod +from .acquire import footprint as footprint_mod from .acquire import mast as mast_mod from .align import diagnostics as align_mod from .drizzle import combine as combine_mod @@ -64,13 +65,26 @@ def reduce_target( ) downloaded = True # Fully-cached re-runs with references already synced stay offline. + # The jwst pipeline syncs its own references lazily through CRDS_PATH, + # so explicit bestrefs is an HST-observatory step. refs_synced = False - if downloaded or not crds_mod.references_present(cache.references_dir, adapter): + if adapter.observatory == "hst" and ( + downloaded + or not crds_mod.references_present(cache.references_dir, adapter) + ): crds_mod.sync_best_references(exposures) refs_synced = True + # Detector-footprint filter: only exposures covering the target enter + # combination — survey visits span many detectors that never touch it, + # and combining them wastes memory (image3 OOM) and time. + cutout_extent = 0.5 * max(spec.cutout_shape) * spec.final_scale + exposures, skipped = footprint_mod.filter_to_target( + exposures, spec.ra, spec.dec, margin_arcsec=cutout_extent + 15.0 + ) record["acquire"] = { "n_exposures": len(exposures), "exposures": [Path(p).name for p in exposures], + "n_skipped_off_target": len(skipped), "downloaded": downloaded, "references_synced": refs_synced, } @@ -94,45 +108,78 @@ def reduce_target( header = hdul[0].header.copy() wht = fits.getdata(wht_path) exptime = header.get("EXPTIME", header.get("TEXPTIME")) - if exptime is None or exptime <= 0: + # Only the HST noise construction divides by exposure time; the JWST path + # reads propagated ERR and needs no exptime — record it if present, but + # never hard-fail a reduction that doesn't use it. + if adapter.combine_backend != "jwst_image3" and ( + exptime is None or exptime <= 0 + ): raise ValueError(f"mosaic header carries no positive EXPTIME: {exptime}") + exptime = float(exptime) if exptime else 0.0 # -- noise ----------------------------------------------------------------- - noise = rms_mod.noise_map_from( - sci, - wht, - exptime=float(exptime), - correlated_noise_factor=drizzle_prov["correlated_noise_factor"], - ) - record["noise"] = { - "recipe": "R * sqrt(max(sci,0)/exptime + 1/wht)", - "correlated_noise_factor": drizzle_prov["correlated_noise_factor"], - "exptime": float(exptime), - "empirical_background_rms": rms_mod.empirical_background_rms( - sci[np.isfinite(noise)] - ), - } + if adapter.combine_backend == "jwst_image3": + from .noise import jwst_rms as jwst_rms_mod + + err = fits.getdata(drizzle_prov["err_path"]).astype(float) + noise, consistency = jwst_rms_mod.noise_map_from_error( + err, + sci, + correlated_noise_factor=drizzle_prov["correlated_noise_factor"], + ) + record["noise"] = { + "recipe": "R * ERR (propagated by calwebb_image3 resample)", + "correlated_noise_factor": drizzle_prov["correlated_noise_factor"], + "exptime": float(exptime), + **consistency, + } + else: + noise = rms_mod.noise_map_from( + sci, + wht, + exptime=float(exptime), + correlated_noise_factor=drizzle_prov["correlated_noise_factor"], + ) + record["noise"] = { + "recipe": "R * sqrt(max(sci,0)/exptime + 1/wht)", + "correlated_noise_factor": drizzle_prov["correlated_noise_factor"], + "exptime": float(exptime), + "empirical_background_rms": rms_mod.empirical_background_rms( + sci[np.isfinite(noise)] + ), + } # -- psf ------------------------------------------------------------------- from astropy.wcs import WCS 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") + if adapter.observatory == "hst": + # 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") + 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 + # from the level-2 pipeline, so no peak cut is applied. Refinement + # (unit-converted cap) tracked in docs/design/jwst.md open items. + peak_max = None 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 - / max_single_exptime, + peak_max=peak_max, ) psf, psf_full, psf_diag = epsf_mod.build_epsf( sci, stars, spec.psf_shape, spec.psf_full_shape @@ -140,13 +187,34 @@ def reduce_target( record["psf"] = psf_diag # -- package ---------------------------------------------------------------- - data_cut = cutout_mod.cutout_to_fits( - sci, header, spec.ra, spec.dec, spec.cutout_shape, out_dir / "data.fits" + data_cut, data_header, center_xy = cutout_mod.make_cutout( + sci, header, spec.ra, spec.dec, spec.cutout_shape + ) + noise_cut, noise_header, _ = cutout_mod.make_cutout( + noise, header, spec.ra, spec.dec, spec.cutout_shape ) - noise_cut = cutout_mod.cutout_to_fits( - noise, header, spec.ra, spec.dec, spec.cutout_shape, out_dir / "noise_map.fits" + # Isolated dead/rejected pixels are masked-by-noise (recorded); the loud + # failure remains for excessive masking or bad pixels near the lens. + data_cut, noise_cut, mask_diag = rms_mod.mask_isolated_bad_pixels( + data_cut, + noise_cut, + center_xy=center_xy, + pixel_scale=spec.final_scale, + region_name=f"{spec.name} cutout", ) rms_mod.assert_finite_within(noise_cut, f"{spec.name} cutout") + cutout_mod.write_fits(data_cut, data_header, out_dir / "data.fits") + cutout_mod.write_fits(noise_cut, noise_header, out_dir / "noise_map.fits") + record["bad_pixel_policy"] = mask_diag + + # The mosaic-wide WHT uniformity mixes coverage tiers across the full + # union footprint; the science verdict belongs to the cutout region. + from .drizzle.diagnostics import check_weight_uniformity + + wht_cut, _, _ = cutout_mod.make_cutout( + wht, header, spec.ra, spec.dec, spec.cutout_shape + ) + record["drizzle"]["weight_uniformity_cutout"] = check_weight_uniformity(wht_cut) fits.PrimaryHDU(psf.astype(np.float32)).writeto( out_dir / "psf.fits", overwrite=True @@ -158,7 +226,8 @@ def reduce_target( "products": ["data.fits", "noise_map.fits", "psf.fits", "psf_full.fits"], "cutout_shape": list(spec.cutout_shape), "pixel_scale": spec.final_scale, - "data_units": drizzle_prov["drizzle_kwargs"]["final_units"], + # Backend-agnostic: both backends stamp BUNIT on the mosaic header. + "data_units": str(header.get("BUNIT", "unknown")), } provenance_mod.write_reduction_json(out_dir, record) diff --git a/autoreduce/psf/epsf.py b/autoreduce/psf/epsf.py index 662d534..9e066df 100644 --- a/autoreduce/psf/epsf.py +++ b/autoreduce/psf/epsf.py @@ -64,12 +64,20 @@ def build_epsf( size += 1 stars = extract_stars(NDData(sci), positions, size=size) - # extract_stars drops stars whose window overruns the mosaic edge, so the - # minimum-star contract must be re-checked on what actually survived. + # Reject stars whose window contains non-finite pixels (coverage edges, + # DQ holes — routine in JWST mosaics): EPSFBuilder's fitter refuses NaN. + from photutils.psf import EPSFStars + + stars = EPSFStars([s for s in stars.all_stars if np.isfinite(s.data).all()]) + + # extract_stars drops stars whose window overruns the mosaic edge and the + # finite cut above drops more, so the minimum-star contract must be + # re-checked on what actually survived. if len(stars) < MIN_STARS: raise InsufficientStarsError( - f"{len(stars)} stars survived cutout extraction (< {MIN_STARS}); " - f"tier 1 ePSF is not viable — select tier 2 (model PSF) explicitly" + f"{len(stars)} stars survived cutout extraction + finite-window " + f"cut (< {MIN_STARS}); tier 1 ePSF is not viable — select tier 2 " + f"(model PSF) explicitly" ) builder = EPSFBuilder(oversampling=oversampling, maxiters=10, progress_bar=False) diff --git a/autoreduce/psf/stars.py b/autoreduce/psf/stars.py index 078ee45..11d0f24 100644 --- a/autoreduce/psf/stars.py +++ b/autoreduce/psf/stars.py @@ -8,7 +8,7 @@ """ from dataclasses import dataclass -from typing import Tuple +from typing import Optional, Tuple import numpy as np @@ -59,13 +59,20 @@ def find_stars( sci: np.ndarray, selection: StarSelection, target_xy: Tuple[float, float], - peak_max: float, + peak_max: Optional[float], ): - """DAOStarFinder detections filtered through every selection cut.""" + """DAOStarFinder detections filtered through every selection cut. + + ``peak_max=None`` disables the saturation cut (surface-brightness-unit + mosaics where a full-well cap is meaningless; saturated cores arrive + blanked from the upstream pipeline). + """ from astropy.stats import sigma_clipped_stats from photutils.detection import DAOStarFinder - _, median, std = sigma_clipped_stats(sci, sigma=3.0) + # Mosaics carry NaN outside the coverage footprint (JWST especially). + finite = np.isfinite(sci) + _, median, std = sigma_clipped_stats(sci[finite], sigma=3.0) finder = DAOStarFinder( fwhm=selection.fwhm_pix, threshold=selection.detection_sigma * std, @@ -75,7 +82,7 @@ def find_stars( roundhi=selection.round_limit, peakmax=peak_max, ) - sources = finder(sci - median) + sources = finder(np.nan_to_num(sci - median), mask=~finite) if sources is None or len(sources) == 0: return None diff --git a/docs/design/jwst.md b/docs/design/jwst.md new file mode 100644 index 0000000..2711e47 --- /dev/null +++ b/docs/design/jwst.md @@ -0,0 +1,81 @@ +# JWST/NIRCam — per-stage deltas vs the HST design + +Phase 3. The first non-HST observatory, and the phase that forced the +**backend dispatch**: stage 3 now routes through +`InstrumentAdapter.combine_backend` (`astrodrizzle` | `jwst_image3`) and the +CRDS server/env shape is adapter-owned (`observatory`, `crds_server_url`; +the jwst pipeline reads `CRDS_PATH` directly — no `jref`-style variable). + +| Stage | Delta vs HST | +|-------|--------------| +| acquire | level-2 **`_cal`** products (calwebb_image2 output, **MJy/sr**), `obs_collection="JWST"`; no explicit bestrefs — the jwst pipeline syncs references lazily through `CRDS_PATH`/jwst-crds | +| align | tweakreg runs *inside* calwebb_image3 (defaults-first); the standalone stage only records WCS provenance | +| combine | `calwebb_image3` (tweakreg / skymatch / outlier_detection / **resample** — the drizzle analogue). The lensing dials map to the resample step: `pixel_scale`, `pixfrac`, `kernel`, `rotation=0`, `weight_type=ivm`. The multi-extension `_i2d` is normalized to standalone sci/wht/err files so downstream stages stay backend-agnostic | +| noise | **read, don't construct**: the resampled `ERR` array (Poisson + read noise + flat, propagated by the pipeline) × the same Casertano R (resample correlates pixels exactly as drizzle does). A blank-sky consistency check (`sky_over_err_floor`) is recorded; disagreement is investigated, never absorbed | +| units | native **MJy/sr** kept (defaults-first — no conversion unless parity demands one); `BUNIT` rides the cutout header | +| psf | tier-1 ePSF unchanged (NaN-masked star finding); **no full-well peak cut** — meaningless in surface-brightness units, and saturated cores arrive blanked from level 2. STPSF is the designated tier-2 back-end (open item) | +| scales | SW native 0.031″ → recommended 0.03″; LW native 0.063″ → recommended 0.06″ (the COSMOS-Web mosaic convention the parity anchor uses). Filter→channel routing via `nircam_adapter_for_filter` | + +## Parity interpretation (maintainer guidance, 2026-07-09) + +The demo dataset descends from the **bespoke COSMOS-Web team pipeline** +(custom 1/f destriping, wisp/snowball handling, their calibration vintage, +mosaics at 0.03″ SW / 0.06″ LW). The acceptance bar is therefore **"close + +internally consistent," not reproduction** — strong lensing needs its own +pipeline (this one), and order-unity data/noise ratios against the team +products are expected and acceptable. What must hold: our own internal +closures (sky vs ERR floor, WHT uniformity over the cutout, masked-pixel +policy) and cross-band consistency of any global scale offset. + +## Validation anchor — the COSMOS-Web ring, four bands + +The autolens_assistant demo dataset +(`dataset/imaging/cosmos_web_ring/wavebands/{F115W,F150W,F277W,F444W}`) +carries modeling-ready products for the ring (RA 150.10048, +1.89301; +[Mercier et al. 2024](https://arxiv.org/abs/2309.15986)) in all four +COSMOS-Web bands — SW at 0.03″/pix (419²), LW at 0.06″/pix (209²), stripped +headers as usual. `scripts/reduce_cosmos_web_ring.py --band ` reduces each +band from MAST `_cal` exposures and reports sub-pixel-registered data/noise +ratios against the demo products (the SLACS-parity method). + +## PSF options — what the JWST weak-lensing and AGN literature says (2026-07-08) + +- **Weak lensing (COSMOS-Web's own practice):** + [ShOpt.jl](https://arxiv.org/abs/2401.11625) (Berman & McCleary 2024) is + COSMOS-Web's PSF characterization tool, benchmarked against **PSFEx** and + **PIFF** on real + simulated COSMOS-Web NIRCam imaging; all model the PSF + **empirically from field stars with low-order polynomial spatial variation + in (X, Y)** across the resampled mosaic. NIRCam PSFs vary with time, + bandpass and field position, so star-based per-mosaic models are the norm. +- **AGN decomposition:** [Zhuang & Shen + 2024](https://arxiv.org/abs/2304.13776) characterize NIRCam PSFs in 8 + filters: spatial FWHM variation shrinks strongly with wavelength (max/RMS + ~20%/5% at F070W → **~3%/0.6% at F444W**); among SWarp / photutils / PSFEx + they find **PSFEx best**; PSF mismatch biases host fluxes high. COSMOS-Web + AGN work ([Zhuang et al. 2024](https://iopscience.iop.org/article/10.3847/1538-4357/ad1517)) + and the galight PSF-library approach (Ding et al.; SHELLQs-JWST) use + curated star libraries / hybrid empirical PSFs; **pure STPSF (WebbPSF) + models are consistently disfavoured vs empirical** for decomposition work. + +**Adopted tiering for JWST (revision of the HST-era tier 2):** + +| Tier | Method | When | +|------|--------|------| +| 1 | single ePSF from mosaic stars (current photutils implementation) | **LW bands** (F277W/F444W): spatial variation ≲1% RMS — a single ePSF at the lens position is adequate for lens-galaxy work | +| 2 | **spatially-varying empirical model evaluated at the lens position** — PSFEx-style polynomial (PSFEx or ShOpt back-end) | **SW bands** (F115W/F150W: ~5% RMS variation) and any weak-lensing-grade use; photutils ranks below PSFEx in the Zhuang & Shen benchmark, so this is the quality upgrade path | +| 2b | STPSF model PSF | fallback only when the field lacks stars — flagged in provenance, never silent (the literature's consistent verdict: empirical beats model for decomposition) | +| 3 | STARRED / PSFr iterative reconstruction | lensed quasars/AGN, unchanged from the HST design | + +Phase 3 ships tier 1; tiers 2/2b are the follow-up (PSFEx/ShOpt are external +binaries/Julia — an integration decision for a dedicated prompt). + +## Open items + +- Tier-2 spatially-varying PSF back-end (PSFEx or ShOpt — see table above); + STPSF as explicit 2b fallback; unit-aware saturation cut for star + selection in MJy/sr mosaics. +- jwst pinned at **1.14.0** by the PyAuto env constraints (astropy 6.1.2); + provenance records the version — revisit when the env's astropy moves. +- COSMOS-Web official reduction (Franco et al.) applies additional + corrections (1/f striping, wisps, snowballs) beyond default calwebb; + parity ratios will show whether they matter at lens-cutout scale. diff --git a/docs/design/roadmap.md b/docs/design/roadmap.md index deb9a91..48485e8 100644 --- a/docs/design/roadmap.md +++ b/docs/design/roadmap.md @@ -22,16 +22,13 @@ adapter #1; nothing outside `instruments/` may mention a detector by name. - Other ACS/WFC3 filters (F435W, F606W…) are config, not code: the adapter already parameterizes the filter-dependent pieces. -## JWST (NIRCam first) - -- Adapter over the STScI `jwst` pipeline: `_cal` products from stage 2 - (`calwebb_image2`), combination via stage 3 (`calwebb_image3`) whose - `resample` step is the drizzle analogue — same defaults-first principle. -- Noise: JWST products carry `ERR`/`VAR_POISSON`/`VAR_RNOISE` arrays, so - stage 4 becomes a *read + resample-consistency check* rather than a - construction — a good test of the stage abstraction. -- PSF: STPSF (formerly WebbPSF) replaces TinyTim in tier 2; tiers 1 and 3 - carry over unchanged. +## JWST (NIRCam first) — **in progress (phase 3, PyAutoReduce#6)** + +- Design deltas live in [`jwst.md`](jwst.md); adapters `nircam_sw`/`nircam_lw` + + the combine-backend dispatch (`astrodrizzle` | `jwst_image3`) implemented; + noise = read propagated ERR × R; validated on the COSMOS-Web ring, four + bands, against the autolens_assistant demo dataset. +- PSF: STPSF tier-2 back-end still open (tier-1 ePSF carries over). ## Per-exposure frame products (`_flt`/`_flc` with cosmic rays) diff --git a/scripts/reduce_cosmos_web_ring.py b/scripts/reduce_cosmos_web_ring.py new file mode 100644 index 0000000..81220b1 --- /dev/null +++ b/scripts/reduce_cosmos_web_ring.py @@ -0,0 +1,145 @@ +""" +JWST integration + acceptance (issue #6): the COSMOS-Web ring, four bands. + +Reduces the ring (RA 150.10048, +1.89301; Mercier et al. 2024) from MAST +level-2 ``_cal`` exposures through calwebb_image3, then compares data/noise +against the autolens_assistant demo dataset for that band (sub-pixel +registered ratios — the SLACS parity method). SW bands (F115W/F150W) output +0.03"/pix; LW (F277W/F444W) 0.06"/pix, matching the demo convention. + +Run: ~/venv/PyAuto/bin/python scripts/reduce_cosmos_web_ring.py --band F444W +Network + jwst pipeline 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.instruments import nircam_adapter_for_filter # noqa: E402 + +RA, DEC = 150.10048, 1.89301 +CACHE_ROOT = REPO / "scripts" / "cache" +OUTPUT_ROOT = REPO / "scripts" / "output" +DEMO_ROOT = Path( + "/home/jammy/Code/PyAutoLabs/autolens_assistant/dataset/imaging/" + "cosmos_web_ring/wavebands" +) + +BANDS = ("F115W", "F150W", "F277W", "F444W") + + +def spec_for(band: str) -> TargetSpec: + adapter = nircam_adapter_for_filter(band) + # Demo cutouts: SW 419x419 @0.03 (12.57"), LW 209x209 @0.06 (12.54") — + # match the demo shapes exactly so parity is pixel-to-pixel. + shape = (419, 419) if adapter.key == "nircam_sw" else (209, 209) + return TargetSpec( + name=f"cosmos_web_ring_{band.lower()}", + ra=RA, + dec=DEC, + instrument=adapter.key, + filter_name=band, + # COSMOS-Web only: the demo parity products are built from program + # 1727 mosaics; other programs at these coords (e.g. 5893) would + # change the depth and skew the noise parity. + proposal_ids=("1727",), + final_scale=adapter.recommended_final_scale, + final_pixfrac=1.0, # COSMOS-Web mosaics use the full drop + cutout_shape=shape, + ) + + +def subpixel_offset(a, b): + a0 = np.nan_to_num(a - np.nanmedian(a)) + b0 = np.nan_to_num(b - np.nanmedian(b)) + corr = np.fft.fftshift( + np.fft.irfft2(np.fft.rfft2(a0) * np.conj(np.fft.rfft2(b0)), s=a0.shape) + ) + peak = np.unravel_index(np.argmax(corr), corr.shape) + + def parabolic(axis): + prev = list(peak); prev[axis] -= 1 + nxt = list(peak); nxt[axis] += 1 + if min(prev[axis], 0) < 0 or nxt[axis] >= corr.shape[axis]: + return 0.0 + cm, c0, cp = corr[tuple(prev)], corr[peak], corr[tuple(nxt)] + denom = cm - 2 * c0 + cp + return 0.0 if denom == 0 else 0.5 * (cm - cp) / denom + + return ( + peak[0] - a.shape[0] // 2 + parabolic(0), + peak[1] - a.shape[1] // 2 + parabolic(1), + ) + + +def compare(band: str, out_dir: Path) -> dict: + from astropy.io import fits + from scipy.ndimage import shift as nd_shift + + new_data = fits.getdata(out_dir / "data.fits").astype(float) + new_noise = fits.getdata(out_dir / "noise_map.fits").astype(float) + demo_dir = DEMO_ROOT / band + demo_data = fits.getdata(demo_dir / "data.fits").astype(float) + demo_noise = fits.getdata(demo_dir / "noise_map.fits").astype(float) + + if new_data.shape != demo_data.shape: + raise ValueError( + f"{band}: shape mismatch new {new_data.shape} vs demo {demo_data.shape}" + ) + + dy, dx = subpixel_offset(demo_data, new_data) + new_data_r = nd_shift(np.nan_to_num(new_data), (dy, dx), order=3) + new_noise_r = nd_shift(np.nan_to_num(new_noise), (dy, dx), order=1) + + bright = demo_data > 10 * np.nanmedian(demo_noise) + data_ratio = new_data_r[bright] / demo_data[bright] + # Exclude masked-by-noise pixels (1e8) and their shift-interpolation + # bleed from the parity statistics. + valid = new_noise_r < 1.0e6 + noise_ratio = np.where(valid, new_noise_r / demo_noise, np.nan) + return { + "band": band, + "offset": [float(dy), float(dx)], + "n_bright": int(bright.sum()), + "data_ratio_median": float(np.nanmedian(data_ratio)), + "noise_ratio_median": float(np.nanmedian(noise_ratio)), + "noise_ratio_16_84": [ + float(np.nanpercentile(noise_ratio, 16)), + float(np.nanpercentile(noise_ratio, 84)), + ], + } + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--band", required=True, choices=[*BANDS, "all"]) + args = parser.parse_args() + bands = BANDS if args.band == "all" else (args.band,) + + for band in bands: + spec = spec_for(band) + record = reduce_target(spec, cache_root=CACHE_ROOT, output_root=OUTPUT_ROOT) + out_dir = OUTPUT_ROOT / spec.name + summary = { + "n_exposures": record["acquire"]["n_exposures"], + "weight_uniformity": record["drizzle"]["weight_uniformity"], + "weight_uniformity_cutout": record["drizzle"].get("weight_uniformity_cutout"), + "correlated_noise_factor": record["noise"]["correlated_noise_factor"], + "sky_over_err_floor": record["noise"].get("sky_over_err_floor"), + "psf": record["psf"], + "parity": compare(band, out_dir), + } + print(f"[{band}] ---- validation ----") + print(json.dumps(summary, indent=2)) + (out_dir / "validation_summary.json").write_text(json.dumps(summary, indent=2)) + + +if __name__ == "__main__": + main() diff --git a/test_autoreduce/test_jwst.py b/test_autoreduce/test_jwst.py new file mode 100644 index 0000000..d469fa6 --- /dev/null +++ b/test_autoreduce/test_jwst.py @@ -0,0 +1,161 @@ +import numpy as np +import pytest + +from autoreduce import instruments +from autoreduce.instruments.nircam import nircam_adapter_for_filter +from autoreduce.noise.jwst_rms import noise_map_from_error + + +class TestNIRCamAdapters: + def test_channels_registered(self): + assert "nircam_sw" in instruments.registered_keys() + assert "nircam_lw" in instruments.registered_keys() + + def test_jwst_routing_fields(self): + for key in ("nircam_sw", "nircam_lw"): + a = instruments.get(key) + assert a.observatory == "jwst" + assert a.combine_backend == "jwst_image3" + assert a.crds_server_url == "https://jwst-crds.stsci.edu" + assert a.mast_obs_collection == "JWST" + assert a.calibrated_suffix == "CAL" + assert not a.supports_cte_correction + + def test_hst_adapters_unchanged(self): + # Phase-1/2 regression: HST adapters keep the astrodrizzle backend + # and the HST CRDS server via the new defaulted fields. + for key in ("acs_wfc", "wfc3_uvis", "wfc3_ir"): + a = instruments.get(key) + assert a.observatory == "hst" + assert a.combine_backend == "astrodrizzle" + assert a.crds_server_url == "https://hst-crds.stsci.edu" + assert a.mast_obs_collection == "HST" + + def test_scales_match_cosmos_web_convention(self): + assert instruments.get("nircam_sw").recommended_final_scale == 0.03 + assert instruments.get("nircam_lw").recommended_final_scale == 0.06 + + def test_filter_routing(self): + assert nircam_adapter_for_filter("F115W").key == "nircam_sw" + assert nircam_adapter_for_filter("F150W").key == "nircam_sw" + assert nircam_adapter_for_filter("F277W").key == "nircam_lw" + assert nircam_adapter_for_filter("f444w").key == "nircam_lw" + with pytest.raises(KeyError, match="unknown NIRCam filter"): + nircam_adapter_for_filter("F814W") + + +class TestBackendDispatch: + def test_unknown_backend_rejected(self): + from dataclasses import replace + + from autoreduce.drizzle.combine import combine + from autoreduce.target import TargetSpec + + broken = replace(instruments.get("acs_wfc"), combine_backend="magic") + spec = TargetSpec(name="x", ra=0.0, dec=0.0) + with pytest.raises(ValueError, match="unknown combine backend"): + combine([], spec, broken, output_dir="/tmp/nowhere") + + +class TestJWSTCRDSEnvironment: + def test_jwst_env_uses_jwst_server_and_no_iraf_var(self, tmp_path, monkeypatch): + import os + + from autoreduce.acquire.crds import configure_environment + + monkeypatch.delenv("CRDS_SERVER_URL", raising=False) + monkeypatch.delenv("CRDS_PATH", raising=False) + env = configure_environment(tmp_path, instruments.get("nircam_lw")) + assert env["CRDS_SERVER_URL"] == "https://jwst-crds.stsci.edu" + assert env["CRDS_PATH"] == str(tmp_path) + # No jref/iref-style variable for the jwst pipeline. + assert set(env) == {"CRDS_SERVER_URL", "CRDS_PATH"} + assert os.environ["CRDS_SERVER_URL"] == "https://jwst-crds.stsci.edu" + + def test_hst_env_still_sets_iraf_var(self, tmp_path, monkeypatch): + from autoreduce.acquire.crds import configure_environment + + monkeypatch.delenv("jref", raising=False) + env = configure_environment(tmp_path, instruments.get("acs_wfc")) + assert env["CRDS_SERVER_URL"] == "https://hst-crds.stsci.edu" + assert env["jref"].endswith("references/hst/acs/") + + +class TestJWSTNoise: + def test_reads_err_and_applies_r(self): + err = np.full((10, 10), 0.02) + sci = np.random.default_rng(0).normal(0.0, 0.02, (10, 10)) + noise, consistency = noise_map_from_error(err, sci, correlated_noise_factor=1.5) + assert noise == pytest.approx(np.full((10, 10), 0.03)) + assert consistency["err_5th_percentile_pre_R"] == pytest.approx(0.02) + + def test_bad_err_pixels_become_nan(self): + err = np.array([[0.02, 0.0], [np.nan, 0.02]]) + sci = np.zeros((2, 2)) + noise, _ = noise_map_from_error(err, sci) + assert np.isnan(noise[0, 1]) and np.isnan(noise[1, 0]) + + def test_shape_mismatch_raises(self): + with pytest.raises(ValueError, match="shape mismatch"): + noise_map_from_error(np.zeros((2, 2)), np.zeros((3, 3))) + + def test_sub_unity_r_raises(self): + with pytest.raises(ValueError): + noise_map_from_error( + np.ones((2, 2)), np.ones((2, 2)), correlated_noise_factor=0.5 + ) + + +def test_find_stars_handles_nan_borders(): + from autoreduce.psf.stars import StarSelection, find_stars + + rng = np.random.default_rng(2) + sci = rng.normal(0.0, 0.01, (200, 200)) + sci[:20, :] = np.nan # coverage border + yy, xx = np.mgrid[0:200, 0:200] + for x0, y0 in [(60, 60), (140, 150), (100, 60)]: + sci += 5.0 * np.exp(-(((xx - x0) ** 2 + (yy - y0) ** 2) / (2 * 1.2**2))) + stars = find_stars( + sci, + StarSelection(min_separation_pix=10.0, edge_margin_pix=25), + target_xy=(0.0, 0.0), + peak_max=None, + ) + assert stars is not None and len(stars) == 3 + + +class TestFootprintFilter: + def _cal_file(self, tmp_path, name, crval): + import numpy as np + from astropy.io import fits + from astropy.wcs import WCS + + wcs = WCS(naxis=2) + wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"] + wcs.wcs.crval = list(crval) + wcs.wcs.crpix = [50.5, 50.5] + wcs.wcs.cdelt = [-1.0 / 3600.0, 1.0 / 3600.0] # 100" x 100" footprint + sci = fits.ImageHDU(np.zeros((100, 100), dtype="f4"), header=wcs.to_header()) + sci.name = "SCI" + path = tmp_path / name + fits.HDUList([fits.PrimaryHDU(), sci]).writeto(path) + return path + + def test_keeps_covering_drops_off_target(self, tmp_path): + from autoreduce.acquire.footprint import filter_to_target + + on = self._cal_file(tmp_path, "on_cal.fits", (150.1, 1.893)) + off = self._cal_file(tmp_path, "off_cal.fits", (150.5, 2.4)) + covering, skipped = filter_to_target( + [on, off], ra=150.1, dec=1.893, margin_arcsec=10.0 + ) + assert covering == [on] and skipped == [off] + + def test_nothing_covering_is_loud(self, tmp_path): + import pytest as _pytest + + from autoreduce.acquire.footprint import filter_to_target + + off = self._cal_file(tmp_path, "off_cal.fits", (150.5, 2.4)) + with _pytest.raises(LookupError, match="cover"): + filter_to_target([off], ra=150.1, dec=1.893, margin_arcsec=10.0) diff --git a/test_autoreduce/test_noise.py b/test_autoreduce/test_noise.py index ba0c9b2..3d02c38 100644 --- a/test_autoreduce/test_noise.py +++ b/test_autoreduce/test_noise.py @@ -109,3 +109,62 @@ def test_empirical_background_rms_recovers_sigma(): rng = np.random.default_rng(0) sci = rng.normal(0.0, 0.02, size=(200, 200)) assert empirical_background_rms(sci) == pytest.approx(0.02, rel=0.05) + + +class TestBadPixelPolicy: + def _pair(self): + data = np.ones((100, 100)) + noise = np.full((100, 100), 0.01) + return data, noise + + def test_clean_cutout_untouched(self): + from autoreduce.noise.rms import mask_isolated_bad_pixels + + data, noise = self._pair() + d, n, diag = mask_isolated_bad_pixels(data, noise, (50.0, 50.0), 0.06) + assert diag["n_masked_pixels"] == 0 + assert (n == 0.01).all() + + def test_isolated_far_pixel_masked_and_recorded(self): + from autoreduce.noise.rms import MASKED_NOISE_VALUE, mask_isolated_bad_pixels + + data, noise = self._pair() + noise[5, 5] = np.nan + d, n, diag = mask_isolated_bad_pixels(data, noise, (50.0, 50.0), 0.06) + assert diag["n_masked_pixels"] == 1 + assert n[5, 5] == MASKED_NOISE_VALUE and d[5, 5] == 0.0 + assert n[0, 0] == 0.01 # rest untouched + + def test_too_many_bad_pixels_is_loud(self): + from autoreduce.noise.rms import mask_isolated_bad_pixels + + data, noise = self._pair() + # 100 isolated singletons on a grid: 1% > 0.5%, no clustering. + noise[::10, ::10] = 0.0 + with pytest.raises(ValueError, match="policy limit"): + mask_isolated_bad_pixels(data, noise, (50.0, 50.0), 0.06) + + def test_bad_pixel_near_lens_is_loud(self): + from autoreduce.noise.rms import mask_isolated_bad_pixels + + data, noise = self._pair() + noise[51, 52] = np.nan # ~0.13" from centre at 0.06"/pix + with pytest.raises(ValueError, match="lens region"): + mask_isolated_bad_pixels(data, noise, (50.0, 50.0), 0.06) + + def test_structured_cluster_is_loud_even_when_small(self): + from autoreduce.noise.rms import mask_isolated_bad_pixels + + data, noise = self._pair() + noise[10:13, 10:13] = np.nan # 3x3 blob: 9 px = 0.09%, but structured + with pytest.raises(ValueError, match="structured"): + mask_isolated_bad_pixels(data, noise, (50.0, 50.0), 0.06) + + def test_scattered_pair_still_maskable(self): + from autoreduce.noise.rms import mask_isolated_bad_pixels + + data, noise = self._pair() + noise[10, 10] = np.nan + noise[10, 11] = np.nan # a pair: each has one bad neighbour + d, n, diag = mask_isolated_bad_pixels(data, noise, (50.0, 50.0), 0.06) + assert diag["n_masked_pixels"] == 2