feat: JWST/NIRCam via the COSMOS-Web ring, four bands (phase 3)

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
 

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

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",

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
 

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

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

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."""

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")

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