diff --git a/CHANGES.rst b/CHANGES.rst index 6fbf175980..e19188168f 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -49,6 +49,12 @@ associations - Exclude NIRISS SOSS data taken with uncalibrated filter F277W from spec2 and tso3 associations. [#8549] +badpix_selfcal +-------------- + +- Added new optional step ``badpix_selfcal`` to the spec2 pipeline to self-calibrate + bad pixels in IFU data. [#8500] + combine_1d ---------- @@ -204,6 +210,9 @@ pipeline virtual ("v") slits and the construction of output file names for each type. [#8442] +- Added new optional step ``badpix_selfcal`` to the ``calwebb_spec2`` to self-calibrate + bad pixels in IFU data. [#8500] + pixel_replace ------------- diff --git a/docs/jwst/badpix_selfcal/arguments.rst b/docs/jwst/badpix_selfcal/arguments.rst new file mode 100644 index 0000000000..98ee18d54b --- /dev/null +++ b/docs/jwst/badpix_selfcal/arguments.rst @@ -0,0 +1,18 @@ +Step Arguments +============== +The ``badpix_selfcal`` step has the following optional arguments. + +``--flagfrac_lower`` (float, default=0.001) + The fraction of pixels to flag as outliers on the low-flux + side of the smoothed-subtracted image. + +``--flagfrac_upper`` (float, default=0.001) + The fraction of pixels to flag as outliers on the high-flux + side of the smoothed-subtracted image. + +``--kernel_size`` (integer, default=15) + The size of the kernel to use for the median filter, which is applied + in the spectral direction to make the smoothed image. + +``--save_flagged_bkg`` (boolean, default=False) + Whether to save the flagged background exposures to fits files. diff --git a/docs/jwst/badpix_selfcal/description.rst b/docs/jwst/badpix_selfcal/description.rst new file mode 100644 index 0000000000..fe68e21e64 --- /dev/null +++ b/docs/jwst/badpix_selfcal/description.rst @@ -0,0 +1,60 @@ +Description +=========== + +:Class: `jwst.badpix_selfcal.BadpixSelfcalStep` +:Alias: badpix_selfcal + +Overview +-------- +The ``badpix_selfcal`` step flags bad pixels in the input data using a self-calibration +technique based on median filtering along the spectral axis. +When additional exposures are available, those are used in combination with the science +exposure to identify bad pixels; when unavailable, the step will be skipped with a warning +unless the ``force_single`` parameter is set True. In that case, the science data alone is +used as its own "background". +This correction is applied to `~jwst.datamodels.IFUImageModel` data +directly after the :ref:`assign_wcs ` correction has been applied +in the :ref:`calwebb_spec2 ` pipeline. + +Input details +------------- +The input data must be in the form of a `~jwst.datamodels.IFUImageModel` or +a `~jwst.datamodels.ModelContainer` containing exactly one +science exposure and any number of additional exposures. +A fits or association file +that can be read into one of these data models is also acceptable. +Any exposure with the metadata attribute ``asn.exptype`` set to +``background`` or ``selfcal`` will be used in conjunction with the science +exposure to construct the combined background image. + +Algorithm +--------- +The algorithm relies on the assumption that bad pixels are outliers in the data along +the spectral axis. The algorithm proceeds as follows: + +* A combined background image is created. If additional (``selfcal`` or ``background``) + exposures are available, + the pixelwise minimum of all background, selfcal, and science exposures is taken. + If no additional exposures are available, the science data itself is passed in + without modification, serving as the "background image" for the rest of the procedure, + i.e., true self-calibration. +* The combined background image is median-filtered, ignoring NaNs, along the spectral axis + with a user-specified kernel size. The default kernel size is 15 pixels. +* The difference between the original background image and the median-filtered background image + is taken. The highest- and lowest-flux pixels in this difference image are + flagged as bad pixels. The default fraction of pixels to flag is 0.1% of the total number of pixels + on each of the high-flux and low-flux ends of the distribution. These fractions can be adjusted + using the ``flagfrac_lower`` and ``flagfrac_upper`` parameters for the low- and high-flux ends + of the distribution, respectively. The total fraction of flagged pixels is thus + ``flagfrac_lower + flagfrac_upper``. +* The bad pixels are flagged in the input data by setting the DQ flag to + "OTHER_BAD_PIXEL" and "DO_NOT_USE". +* The bad pixels are also flagged in each exposure with ``asn.exptype`` equal to ``background``, + if available. + +Output product +-------------- +The output is a new copy of the input `~jwst.datamodels.IFUImageModel`, with the +bad pixels flagged. If the entire ``calwebb_spec2`` pipeline is run, the background +exposures passed into the ``background`` step will include the flags from the +``badpix_selfcal`` step. diff --git a/docs/jwst/badpix_selfcal/index.rst b/docs/jwst/badpix_selfcal/index.rst new file mode 100644 index 0000000000..74fc4d2296 --- /dev/null +++ b/docs/jwst/badpix_selfcal/index.rst @@ -0,0 +1,14 @@ +.. _badpix_selfcal_step: + +========================== +Bad Pixel Self-Calibration +========================== + +.. toctree:: + :maxdepth: 2 + + description.rst + arguments.rst + reference_files.rst + +.. automodapi:: jwst.badpix_selfcal diff --git a/docs/jwst/badpix_selfcal/reference_files.rst b/docs/jwst/badpix_selfcal/reference_files.rst new file mode 100644 index 0000000000..b433d1288f --- /dev/null +++ b/docs/jwst/badpix_selfcal/reference_files.rst @@ -0,0 +1,4 @@ +Reference Files +=============== +The ``badpix_selfcal`` step does not use any reference files. + diff --git a/docs/jwst/package_index.rst b/docs/jwst/package_index.rst index 9a75211246..c5d2096ae8 100644 --- a/docs/jwst/package_index.rst +++ b/docs/jwst/package_index.rst @@ -13,6 +13,7 @@ Package Index associations/index.rst background_step/index.rst background_subtraction/index.rst + badpix_selfcal/index.rst barshadow/index.rst charge_migration/index.rst combine_1d/index.rst diff --git a/docs/jwst/pipeline/calwebb_spec2.rst b/docs/jwst/pipeline/calwebb_spec2.rst index fa53774c89..af3284c6c6 100644 --- a/docs/jwst/pipeline/calwebb_spec2.rst +++ b/docs/jwst/pipeline/calwebb_spec2.rst @@ -47,6 +47,8 @@ TSO exposures. The instrument mode abbreviations used in the table are as follow +==========================================================+=====+=====+=====+=====+=====+=====+=================+======+========+=====+ | :ref:`assign_wcs ` | |c| | |c| | |c| | |c| | |c| | |c| | |c| | |c| | |c| | |c| | +----------------------------------------------------------+-----+-----+-----+-----+-----+-----+-----------------+------+--------+-----+ +| :ref:`badpix_selfcal `\ :sup:`2` | | | |c| | | | |c| | | | | | ++----------------------------------------------------------+-----+-----+-----+-----+-----+-----+-----------------+------+--------+-----+ | :ref:`msaflagopen ` | | |c| | |c| | | | | | | | | +----------------------------------------------------------+-----+-----+-----+-----+-----+-----+-----------------+------+--------+-----+ | :ref:`nsclean ` | |c| | |c| | |c| | | | | | | | | diff --git a/jwst/badpix_selfcal/__init__.py b/jwst/badpix_selfcal/__init__.py new file mode 100644 index 0000000000..2df875d38b --- /dev/null +++ b/jwst/badpix_selfcal/__init__.py @@ -0,0 +1,4 @@ +from .badpix_selfcal_step import BadpixSelfcalStep + + +__all__ = ['BadpixSelfcalStep'] \ No newline at end of file diff --git a/jwst/badpix_selfcal/badpix_selfcal.py b/jwst/badpix_selfcal/badpix_selfcal.py new file mode 100644 index 0000000000..f6724a8465 --- /dev/null +++ b/jwst/badpix_selfcal/badpix_selfcal.py @@ -0,0 +1,91 @@ +from __future__ import annotations +import numpy as np +import jwst.datamodels as dm +from jwst.outlier_detection.outlier_detection_ifu import medfilt +from stdatamodels.jwst.datamodels.dqflags import pixel +import warnings + + +def badpix_selfcal(minimg: np.ndarray, + flagfrac_lower: float = 0.001, + flagfrac_upper: float = 0.001, + kernel_size: int = 15, + dispaxis=None) -> np.ndarray: + """ + Flag residual artifacts as bad pixels in the DQ array of a JWST exposure + + Parameters + ---------- + minimg : np.ndarray + Selfcal data of shape (x, y), i.e., after some operation has + already been taken to combine multiple exposures, + typically a MIN operation. + flagfrac_lower : float + Fraction of pixels to flag on the low end + flagfrac_upper : float + Fraction of pixels to flag on the high end + kernel_size : int + Size of kernel for median filter + dispaxis : int + Dispersion axis, either 1 or 2. If None, a two-dimensional + median filter is applied. + + Returns + ------- + flagged_indices : np.ndarray + Indices of the flagged pixels, + shaped like output from np.where + """ + if dispaxis not in [1, 2, None]: + raise ValueError("dispaxis must be either 1 or 2, or None.") + + # Determine outliers using median filter + if dispaxis is None: + kern_size = (kernel_size, kernel_size) + elif dispaxis == 2: + # check that this is the right way around! + kern_size = (kernel_size, 1) + elif dispaxis == 1: + kern_size = (1, kernel_size) + + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", category=RuntimeWarning, message="All-NaN") + smoothed = medfilt(minimg, kern_size) + minimg_hpf = minimg - smoothed + + # Flag outliers using percentile cutoff + flag_low, flag_high = np.nanpercentile(minimg_hpf, [flagfrac_lower * 100, (1 - flagfrac_upper) * 100]) + bad = (minimg_hpf > flag_high) | (minimg_hpf < flag_low) + flagged_indices = np.where(bad) + + return flagged_indices + + +def apply_flags(input_model: dm.IFUImageModel, flagged_indices: np.ndarray) -> dm.IFUImageModel: + """ + Apply the flagged indices to the input model. Sets the flagged pixels to NaN + and the DQ flag to DO_NOT_USE + OTHER_BAD_PIXEL + + Parameters + ---------- + input_model : dm.IFUImageModel + Input science data to be corrected + flagged_indices : np.ndarray + Indices of the flagged pixels, + shaped like output from np.where + + Returns + ------- + output_model : dm.IFUImageModel + Flagged data model + """ + + input_model.dq[flagged_indices] |= pixel["DO_NOT_USE"] + pixel["OTHER_BAD_PIXEL"] + + input_model.data[flagged_indices] = np.nan + input_model.err[flagged_indices] = np.nan + input_model.var_poisson[flagged_indices] = np.nan + input_model.var_rnoise[flagged_indices] = np.nan + input_model.var_flat[flagged_indices] = np.nan + + return input_model diff --git a/jwst/badpix_selfcal/badpix_selfcal_step.py b/jwst/badpix_selfcal/badpix_selfcal_step.py new file mode 100644 index 0000000000..a8cde5d888 --- /dev/null +++ b/jwst/badpix_selfcal/badpix_selfcal_step.py @@ -0,0 +1,209 @@ +import warnings +from ..stpipe import Step +from . import badpix_selfcal +import numpy as np +from jwst import datamodels as dm + +__all__ = ["BadpixSelfcalStep"] + + +class BadpixSelfcalStep(Step): + """ + BadpixSelfcalStep: Flags residual artifacts as bad pixels in the DQ array + of a JWST exposure using a median filter and percentile cutoffs. + + All input exposures in the association file (or manually-provided bkg_list) are combined + into a single background model using a MIN operation. The bad pixels are then identified + using a median filter and percentile cutoffs, and applied to the science data by setting + the flagged pixels, errors, and variances to NaN + and the DQ flag to DO_NOT_USE + OTHER_BAD_PIXEL. + """ + + class_alias = "badpix_selfcal" + + spec = """ + flagfrac_lower = float(default=0.001, min=0.0, max=0.5) # fraction of pixels to flag on the low-flux end + flagfrac_upper = float(default=0.001, min=0.0, max=0.5) # fraction of pixels to flag on the high-flux end + kernel_size = integer(default=15, min=1) # size of kernel for median filter + force_single = boolean(default=False) # force single input exposure + save_flagged_bkg = boolean(default=False) # save flagged background exposures to file + skip = boolean(default=True) + """ + + def save_model(self, model, *args, **kwargs): + """Override save_model to suppress index 0 when save_model is True + """ + kwargs["idx"] = None + return Step.save_model(self, model, *args, **kwargs) + + def save_bkg(self, bkg_list, suffix="badpix_selfcal_bkg"): + """Save the background exposures to file with correct indexing + """ + for i, bkg_model in enumerate(bkg_list): + self.save_model(bkg_model, suffix=suffix + f"_{str(i)}") + + def process(self, input, selfcal_list=None, bkg_list=None): + """ + Flag residual artifacts as bad pixels in the DQ array of a JWST exposure + + Parameters + ---------- + input : JWST data model or association + input science data to be corrected, or tuple of (sci, bkg, selfcal) + + selfcal_list : list of ImageModels or filenames, default None + exposures to include as part of median background model used to find bad pixels, + but that are not flagged and returned as background exposures + + bkg_list : list of ImageModels or filenames, default None + exposures to include as part of median background model used to find bad pixels, + and that are flagged and returned as background exposures + + Returns + ------- + output : JWST data model or association + data model with CRs flagged + + Notes + ----- + If an association file is read in, all exposures in the + association file, including science, background, and selfcal exposures, + are included in the MIN frame from which outliers are detected. + If selfcal_list and/or bkg_list are specified manually, they are appended to any + selfcal or background exposures found in the input association file. + If selfcal_list and bkg_list are both set to None and input is + a single science exposure, the step will be skipped with a warning unless + the force_single parameter is set True. + In that case, the input exposure will be used as the sole background exposure, + i.e., true self-calibration. + """ + input_sci, selfcal_list, bkg_list = _parse_inputs(input, selfcal_list, bkg_list) + + # ensure that there are background exposures to use, otherwise skip the step + # unless forced + if (len(selfcal_list + bkg_list) == 0) and (not self.force_single): + self.log.warning("No selfcal or background exposures provided for self-calibration. " + "Skipping step.") + self.log.info("If you want to force self-calibration with the science " + "exposure alone (generally not recommended), set force_single=True.") + input_sci.meta.cal_step.badpix_selfcal = 'SKIPPED' + return input_sci, bkg_list + + # get the dispersion axis + try: + dispaxis = input_sci.meta.wcsinfo.dispersion_direction + except AttributeError: + self.log.warning("Dispersion axis not found in input science image metadata. " + "Kernel for self-calibration will be two-dimensional.") + dispaxis = None + + # collapse all selfcal exposures into a single background model + # note that selfcal_list includes the science exposure. This is expected. + # all exposures are combined into a single background model using a MIN operation. + selfcal_list = [input_sci] + selfcal_list + selfcal_3d = [] + for i, selfcal_model in enumerate(selfcal_list): + selfcal_3d.append(selfcal_model.data) + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", category=RuntimeWarning, message="All-NaN") + minimg = np.nanmin(np.asarray(selfcal_3d), axis=0) + bad_indices = badpix_selfcal.badpix_selfcal(minimg, self.flagfrac_lower, self.flagfrac_upper, self.kernel_size, dispaxis) + + # apply the flags to the science data + input_sci = badpix_selfcal.apply_flags(input_sci, bad_indices) + + # apply the flags to the background data to be passed to background sub step + if len(bkg_list) > 0: + for i, background_model in enumerate(bkg_list): + bkg_list[i] = badpix_selfcal.apply_flags(dm.open(background_model), bad_indices) + + if self.save_flagged_bkg: + self.save_bkg(bkg_list) + + input_sci.meta.cal_step.badpix_selfcal = 'COMPLETE' + return input_sci, bkg_list + + +def _parse_inputs(input, selfcal_list, bkg_list): + """ + Parse the input to the step. This is a helper function that is used in the + command line interface to the step. + + Parameters + ---------- + input : str + input science exposure or association + + selfcal_list : list of ImageModels or filenames, default None + exposures to include as part of median background model used to find bad pixels, + but that are not flagged and returned as background exposures + + bkg_list : list of ImageModels or filenames, default None + exposures to include as part of median background model used to find bad pixels, + and that are flagged and returned as background exposures + + Returns + ------- + input : JWST data model or association + input science data to be corrected + + selfcal_list : list of ImageModels or filenames to use for selfcal + """ + if selfcal_list is None: + selfcal_list = [] + selfcal_list = [dm.open(selfcal) for selfcal in selfcal_list] + if bkg_list is None: + bkg_list = [] + bkg_list = [dm.open(bkg) for bkg in bkg_list] + selfcal_list = selfcal_list + bkg_list + + with dm.open(input) as input_data: + + # find science and background exposures in association file + if isinstance(input_data, dm.ModelContainer): + + sci_models, bkg_list_asn, selfcal_list_asn = split_container_by_asn_exptype( + input_data, exptypes=['science', 'background', 'selfcal']) + + selfcal_list = selfcal_list + list(bkg_list_asn) + list(selfcal_list_asn) + bkg_list += list(bkg_list_asn) + + if len(sci_models) > 1: + raise ValueError("Input data contains multiple science exposures. " + "This is not supported in calwebb_spec2 steps.") + input_sci = sci_models[0] + + elif isinstance(input_data, dm.IFUImageModel) or isinstance(input_data, dm.ImageModel): + + input_sci = input_data + + else: + raise TypeError("Input data is not a ModelContainer, ImageModel, or IFUImageModel. " + "Cannot continue.") + + return input_sci, selfcal_list, bkg_list + + +def split_container_by_asn_exptype(container: dm.ModelContainer, exptypes: list) -> list: + """ + Split a ModelContainer into a list of ImageModels based on exposure type. + + Parameters + ---------- + container : ModelContainer + input data + + exptype : str + exposure type to split on + + Returns + ------- + split_list : list of lists + lists of ImageModels + """ + split_list = [] + for exptype in exptypes: + good_models = [container[j] for j in range(len(container)) if container[j].meta.asn.exptype == exptype] + split_list.append(good_models) + + return split_list diff --git a/jwst/badpix_selfcal/tests/__init__.py b/jwst/badpix_selfcal/tests/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/jwst/badpix_selfcal/tests/test_badpix_selfcal.py b/jwst/badpix_selfcal/tests/test_badpix_selfcal.py new file mode 100644 index 0000000000..308972ae58 --- /dev/null +++ b/jwst/badpix_selfcal/tests/test_badpix_selfcal.py @@ -0,0 +1,317 @@ +import json +import pytest +import numpy as np +from jwst.badpix_selfcal.badpix_selfcal_step import BadpixSelfcalStep, _parse_inputs +from jwst.badpix_selfcal.badpix_selfcal import badpix_selfcal, apply_flags +from jwst import datamodels as dm +from jwst.assign_wcs import AssignWcsStep, miri +from gwcs import wcs +from astropy.io import fits +from stdatamodels.jwst.datamodels.dqflags import pixel + +wcs_kw = {'wcsaxes': 3, 'ra_ref': 165, 'dec_ref': 54, + 'v2_ref': -8.3942412, 'v3_ref': -5.3123744, 'roll_ref': 37, + 'crpix1': 1024, 'crpix2': 1024, 'crpix3': 0, + 'cdelt1': .08, 'cdelt2': .08, 'cdelt3': 1, + 'ctype1': 'RA---TAN', 'ctype2': 'DEC--TAN', 'ctype3': 'WAVE', + 'pc1_1': 1, 'pc1_2': 0, 'pc1_3': 0, + 'pc2_1': 0, 'pc2_2': 1, 'pc2_3': 0, + 'pc3_1': 0, 'pc3_2': 0, 'pc3_3': 1, + 'cunit1': 'deg', 'cunit2': 'deg', 'cunit3': 'um', + } + +# these are often warm pixels, so need to be much larger than noise, +# which has std dev of 1, +# but smaller than smoothly-varying background, which maxes around 20 +hotpixel_intensity = 10 +outlier_indices = [(100, 100), (300, 300), (500, 600), (1000, 900)] + + +def create_hdul(detector, channel, band): + hdul = fits.HDUList() + phdu = fits.PrimaryHDU() + phdu.header['telescop'] = "JWST" + phdu.header['filename'] = "test" + channel + band + phdu.header['instrume'] = 'MIRI' + phdu.header['detector'] = detector + phdu.header['CHANNEL'] = channel + phdu.header['BAND'] = band + phdu.header['time-obs'] = '8:59:37' + phdu.header['date-obs'] = '2017-09-05' + phdu.header['exp_type'] = 'MIR_MRS' + scihdu = fits.ImageHDU() + scihdu.header['EXTNAME'] = "SCI" + scihdu.header.update(wcs_kw) + hdul.append(phdu) + hdul.append(scihdu) + return hdul + + +def create_reference_files(datamodel): + refs = {} + step = AssignWcsStep() + for reftype in AssignWcsStep.reference_file_types: + refs[reftype] = step.get_reference_file(datamodel, reftype) + + return refs + + +@pytest.fixture(scope="module") +def background(): + """ + Create background IFUImageModel for testing. This is a mockup of the expected + background data in a .rate file. + Three components: random noise, low-order variability, and outliers + """ + + # random noise + rng = np.random.default_rng(seed=77) + shp = (1024, 1032) + noise = rng.standard_normal(shp) + + # make 2-d polynomial representing background level + c = np.array([[1, 3, 5], [2, 4, 6]]) + x = np.linspace(-1, 1, shp[0]) + y = np.linspace(-1, 1, shp[1]) + low_order_variability = np.polynomial.polynomial.polygrid2d(x, y, c) + + # add some outliers + outliers = np.zeros(shp) + for idx in outlier_indices: + outliers[idx] = hotpixel_intensity + # one negative one just for fun + outliers[100, 100] = -hotpixel_intensity + + mock_data = low_order_variability + noise + outliers + + # build an IFUImageModel from these data and give it a wcs + hdul = create_hdul(detector="MIRIFULONG", channel="34", band="LONG") + hdul[1].data = mock_data + + im = dm.IFUImageModel(hdul) + ref = create_reference_files(im) + pipeline = miri.create_pipeline(im, ref) + wcsobj = wcs.WCS(pipeline) + im.meta.wcs = wcsobj + + return im + + +@pytest.fixture(scope="module") +def sci(background): + """Create science IFUImageMoel for testing. + Same data as background with different noise. + """ + hdul = create_hdul(detector="MIRIFULONG", channel="34", band="LONG") + hdul[1].data = background.data.copy() + rng = np.random.default_rng(seed=99) + hdul[1].data += rng.standard_normal(hdul[1].data.shape) + im = dm.IFUImageModel(hdul) + im.meta.wcs = background.meta.wcs + return im + + +@pytest.fixture(scope="module") +def asn(tmp_cwd_module, background, sci): + """Create association for testing. Needs at least + two background images to properly test the step.""" + + sci_path = tmp_cwd_module / "sci.fits" + sci.save(sci_path) + + for root in ["bkg0", "bkg1", "selfcal0", "selfcal1"]: + path = tmp_cwd_module / (root + ".fits") + background.save(path) + + asn_table = { + "asn_pool": "singleton", + "products": [ + { + "members": [ + { + "expname": "sci.fits", + "exptype": "science" + }, + { + "expname": "bkg0.fits", + "exptype": "background" + }, + { + "expname": "bkg1.fits", + "exptype": "background" + }, + { + "expname": "selfcal0.fits", + "exptype": "selfcal" + }, + { + "expname": "selfcal1.fits", + "exptype": "selfcal" + } + ] + } + ] + } + with open(tmp_cwd_module / 'tmp_asn.json', 'w') as f: + json.dump(asn_table, f) + + container = dm.open(tmp_cwd_module / 'tmp_asn.json') + + return container + + +def test_input_parsing(asn, sci, background): + """Test that the input parsing function correctly identifies the science, + background, and selfcal exposures in the association file + given association vs imagemodel inputs, and selfcal_list vs not + + Note that science exposure gets added to selfcal_list later, not in parse_inputs + """ + + # basic association case. Both background and selfcal get into the list + input_sci, selfcal_list, bkg_list = _parse_inputs(asn, [], []) + assert isinstance(input_sci, dm.IFUImageModel) + assert len(bkg_list) == 2 + assert len(selfcal_list) == 4 + + # association with background_list provided + input_sci, selfcal_list, bkg_list = _parse_inputs(asn, [], [background,] * 3) + assert isinstance(input_sci, dm.IFUImageModel) + assert len(bkg_list) == 5 + assert len(selfcal_list) == 7 + + # association with selfcal_list provided + input_sci, selfcal_list, bkg_list = _parse_inputs(asn, [background,] * 3, []) + assert isinstance(input_sci, dm.IFUImageModel) + assert len(bkg_list) == 2 + assert len(selfcal_list) == 7 + + # single science exposure + input_sci, selfcal_list, bkg_list = _parse_inputs(sci, [], []) + assert isinstance(input_sci, dm.IFUImageModel) + assert len(bkg_list) == 0 + assert len(selfcal_list) == 0 + + # single science exposure with selfcal_list and bkg_list provided + input_sci, selfcal_list, bkg_list = _parse_inputs(sci, [background,] * 3, [background,] * 1) + assert isinstance(input_sci, dm.IFUImageModel) + assert len(bkg_list) == 1 + assert len(selfcal_list) == 4 + + +def test_background_flagger_mrs(background): + """ + Ensure the right number of outliers are found, and that + true outliers are among those. + """ + + bg = background.data + + # pass into the MRSBackgroundFlagger and check it found the right pixels + flagfrac = 0.001 + result = badpix_selfcal(bg, flagfrac_lower=flagfrac, flagfrac_upper=flagfrac) + result_tuples = [(i, j) for i, j in zip(*result)] + + # check that the hot pixels were among those flagged + for idx in outlier_indices: + assert idx in result_tuples + + # check that the number of flagged pixels is as expected + assert np.isclose(len(result_tuples) / bg.size, flagfrac * 2, atol=0.0001) + + +def test_apply_flags(background): + """ + Ensure that flagged pixels are set to NaN in the data and err arrays, + and that the DQ flag is set to 1. + """ + + flagged_indices = np.array(outlier_indices).T + flagged_indices = (flagged_indices[0], flagged_indices[1]) + + flagged = apply_flags(background, flagged_indices) + + # check that flagged pixels are NaN in data and err arrays + for idx in outlier_indices: + assert np.isnan(flagged.data[idx]) + assert np.isnan(flagged.err[idx]) + assert np.isnan(flagged.var_poisson[idx]) + assert np.isnan(flagged.var_rnoise[idx]) + assert np.isnan(flagged.var_flat[idx]) + + # check that DQ flag is set properly + for idx in outlier_indices: + assert flagged.dq[idx] == pixel["DO_NOT_USE"] + pixel["OTHER_BAD_PIXEL"] + + +@pytest.mark.parametrize("dset", ["sci", "asn"]) +def test_badpix_selfcal_step(request, dset): + """Test the badpix_selfcal step. This is a functional test that checks + that the step runs without error. The output will be checked by the + regtest. + """ + input_data = request.getfixturevalue(dset) + result = BadpixSelfcalStep.call(input_data, skip=False, force_single=True) + + assert result[0].meta.cal_step.badpix_selfcal == "COMPLETE" + if dset == "sci": + assert len(result) == 2 + assert isinstance(result[0], dm.IFUImageModel) + assert len(result[1]) == 0 + else: + # should return sci, (bkg0, bkg1) but not selfcal0, selfcal1 + assert len(result) == 2 + assert isinstance(result[0], dm.IFUImageModel) + assert len(result[1]) == 2 + + +def test_expected_fail_sci(sci): + """Test that the step fails as expected when given only a science exposure + and force_single is False. + """ + result = BadpixSelfcalStep.call(sci, skip=False, force_single=False) + assert result[0].meta.cal_step.badpix_selfcal == "SKIPPED" + + +@pytest.fixture(scope="module") +def vertical_striping(): + """2-D array with vertical stripes and some outliers, used to test that + the step is flagging in the correct (along-dispersion) direction + """ + shp = (102, 96) + stripes = np.zeros(shp) + stripes[:, ::4] = hotpixel_intensity + + # add some outliers + outliers = np.zeros(shp) + for idx in outlier_indices: + i, j = int(idx[0] / 10), int(idx[1] / 10) + outliers[i, j] = hotpixel_intensity + data = stripes + outliers + + return data + + +@pytest.mark.parametrize("dispaxis", [None, 1, 2]) +def test_dispersion_direction(vertical_striping, dispaxis): + """ + Test that the flagging operates as expected as a function of direction. + + Of the four outliers, only one lies atop a stripe, so the step should only + flag one in four as an outlier unless median filter occurs in the along-stripe + direction. + """ + + flagged_indices = badpix_selfcal(vertical_striping, dispaxis=dispaxis) + + if dispaxis == 2: + + assert len(flagged_indices[0]) == 4 + + elif dispaxis == 1: + + assert len(flagged_indices[0]) == 1 + + elif dispaxis is None: + + assert len(flagged_indices[0]) == 1 diff --git a/jwst/lib/suffix.py b/jwst/lib/suffix.py index cab5be43a9..3fff82ba9f 100644 --- a/jwst/lib/suffix.py +++ b/jwst/lib/suffix.py @@ -195,7 +195,8 @@ 'image2pipeline', 'klip', 'emicorr', - 'emicorrstep' + 'emicorrstep', + 'badpixselfcalstep', } diff --git a/jwst/pipeline/calwebb_spec2.py b/jwst/pipeline/calwebb_spec2.py index aa88fcb87e..57245439b4 100644 --- a/jwst/pipeline/calwebb_spec2.py +++ b/jwst/pipeline/calwebb_spec2.py @@ -13,6 +13,7 @@ # step imports from ..assign_wcs import assign_wcs_step from ..background import background_step +from ..badpix_selfcal import badpix_selfcal_step from ..barshadow import barshadow_step from ..cube_build import cube_build_step from ..extract_1d import extract_1d_step @@ -65,6 +66,7 @@ class Spec2Pipeline(Pipeline): # Define aliases to steps step_defs = { 'assign_wcs': assign_wcs_step.AssignWcsStep, + 'badpix_selfcal': badpix_selfcal_step.BadpixSelfcalStep, 'msa_flagging': msaflagopen_step.MSAFlagOpenStep, 'nsclean': nsclean_step.NSCleanStep, 'bkg_subtract': background_step.BackgroundStep, @@ -244,6 +246,20 @@ def process_exposure_product( # Steps whose order is the same for all types of input: + # Self-calibrate to flag bad/warm pixels, and apply flags + # to both background and science exposures. + # skipped by default for all modes + result = self.badpix_selfcal( + calibrated, members_by_type['selfcal'], members_by_type['background'], + ) + if isinstance(result, datamodels.JwstDataModel): + # if step is skipped, unchanged sci exposure is returned + calibrated = result + else: + # if step actually occurs, then flagged backgrounds are also returned + calibrated, bkg_outlier_flagged = result[0], result[1] + members_by_type['background'] = bkg_outlier_flagged + # apply msa_flagging (flag stuck open shutters for NIRSpec IFU and MOS) calibrated = self.msa_flagging(calibrated) @@ -254,7 +270,6 @@ def process_exposure_product( # If there is only one `imprint` member, this imprint exposure is subtracted from all the # science and background exposures. Otherwise, there will be as many `imprint` members as # there are science plus background members. - calibrated = self.imprint_subtract(calibrated, members_by_type['imprint']) # for each background image subtract an associated leak cal diff --git a/jwst/pipeline/tests/test_calwebb_spec2.py b/jwst/pipeline/tests/test_calwebb_spec2.py index 59b2148bc1..d060aa9cd0 100644 --- a/jwst/pipeline/tests/test_calwebb_spec2.py +++ b/jwst/pipeline/tests/test_calwebb_spec2.py @@ -64,6 +64,7 @@ def run_spec2_pipeline(make_dummy_rate_file, request): and skipping most of the rest for runtime ''' args = ["calwebb_spec2", INPUT_FILE, + "--steps.badpix_selfcal.skip=true", "--steps.msa_flagging.skip=true", "--steps.nsclean.skip=true", "--steps.flat_field.skip=true", @@ -92,6 +93,7 @@ def run_spec2_pipeline_asn(make_dummy_association, request): args = ["calwebb_spec2", INPUT_ASN, f"--logcfg={LOGCFG}", + "--steps.badpix_selfcal.skip=true", "--steps.msa_flagging.skip=true", "--steps.nsclean.skip=true", "--steps.flat_field.skip=true", diff --git a/jwst/regtest/test_miri_mrs_badpix_selfcal.py b/jwst/regtest/test_miri_mrs_badpix_selfcal.py new file mode 100644 index 0000000000..e8b1effb49 --- /dev/null +++ b/jwst/regtest/test_miri_mrs_badpix_selfcal.py @@ -0,0 +1,73 @@ +import pytest +from astropy.io.fits.diff import FITSDiff +from jwst.stpipe import Step +import os + + +@pytest.fixture(scope="module") +def run_pipeline_background(rtdata_module): + + rtdata = rtdata_module + rtdata.get_asn("miri/mrs/jw01204-o021_20240127t024203_spec2_00010_asn.json") + Step.from_cmdline(['calwebb_spec2', rtdata.input, + "--steps.badpix_selfcal.save_results=True", + "--steps.badpix_selfcal.save_flagged_bkg=True", + "--steps.badpix_selfcal.flagfrac_lower=0.0005", + "--steps.badpix_selfcal.skip=False"]) + rtdata.output = "jw01204021001_02101_00004_mirifulong_badpix_selfcal.fits" + return rtdata + + +@pytest.fixture(scope="module") +def run_pipeline_selfcal(rtdata_module): + '''Identical pipeline run to above, but input asn sets all background exposures as `selfcal` type. + ''' + rtdata = rtdata_module + rtdata.get_asn("miri/mrs/jw01204-o021_20240127t024203_spec2_00010_selfcal_asn.json") + Step.from_cmdline(['calwebb_spec2', rtdata.input, + "--steps.badpix_selfcal.save_results=True", + "--steps.badpix_selfcal.flagfrac_lower=0.0005", + "--steps.badpix_selfcal.skip=False"]) + rtdata.output = "jw01204021001_02101_00004_mirifulong_badpix_selfcal.fits" + + return rtdata + + +@pytest.mark.bigdata +def test_miri_mrs_badpix_selfcal(run_pipeline_selfcal, fitsdiff_default_kwargs): + """Run a test for MIRI MRS data with dedicated background exposures.""" + + rtdata = run_pipeline_selfcal + + # Get the truth file + rtdata.get_truth("truth/test_miri_mrs_badpix_selfcal/jw01204021001_02101_00004_mirifulong_badpix_selfcal.fits") + + # Compare the results + diff = FITSDiff(rtdata.output, rtdata.truth, **fitsdiff_default_kwargs) + assert diff.identical, diff.report() + + # check the bkg files in the background case, but not in the selfcal case + for idx in range(4): + fname = f"jw01204021001_02101_00004_mirifulong_badpix_selfcal_bkg_{idx}.fits" + assert not os.path.isfile(fname) + + +@pytest.mark.bigdata +def test_miri_mrs_badpix_selfcal_bkg(run_pipeline_background, fitsdiff_default_kwargs): + """Run a test for MIRI MRS data with dedicated background exposures.""" + + rtdata = run_pipeline_background + + # Get the truth file + rtdata.get_truth("truth/test_miri_mrs_badpix_selfcal/jw01204021001_02101_00004_mirifulong_badpix_selfcal.fits") + + # Compare the results + diff = FITSDiff(rtdata.output, rtdata.truth, **fitsdiff_default_kwargs) + assert diff.identical, diff.report() + + # check the bkg files in the background case, but not in the selfcal case + for idx in range(4): + fname = f"jw01204021001_02101_00004_mirifulong_badpix_selfcal_bkg_{idx}.fits" + truth = rtdata.get_truth(f"truth/test_miri_mrs_badpix_selfcal/{fname}") + diff = FITSDiff(fname, truth, **fitsdiff_default_kwargs) + assert diff.identical, diff.report() \ No newline at end of file diff --git a/jwst/step.py b/jwst/step.py index 7a5d124b59..1c9b67b58f 100644 --- a/jwst/step.py +++ b/jwst/step.py @@ -4,6 +4,7 @@ from .assign_mtwcs.assign_mtwcs_step import AssignMTWcsStep from .assign_wcs.assign_wcs_step import AssignWcsStep from .background.background_step import BackgroundStep +from .badpix_selfcal.badpix_selfcal_step import BadpixSelfcalStep from .barshadow.barshadow_step import BarShadowStep from .charge_migration.charge_migration_step import ChargeMigrationStep from .combine_1d.combine_1d_step import Combine1dStep @@ -68,6 +69,7 @@ "AssignMTWcsStep", "AssignWcsStep", "BackgroundStep", + "BadpixSelfcalStep", "BarShadowStep", "Combine1dStep", "StackRefsStep", diff --git a/jwst/stpipe/integration.py b/jwst/stpipe/integration.py index b43f4abff5..597a6afe3e 100644 --- a/jwst/stpipe/integration.py +++ b/jwst/stpipe/integration.py @@ -37,6 +37,7 @@ def get_steps(): ("jwst.step.AssignMTWcsStep", 'assign_mtwcs', False), ("jwst.step.AssignWcsStep", 'assign_wcs', False), ("jwst.step.BackgroundStep", 'background', False), + ("jwst.step.BadpixSelfcalStep", 'badpix_selfcal', False), ("jwst.step.BarShadowStep", 'barshadow', False), ("jwst.step.Combine1dStep", 'combine_1d', False), ("jwst.step.StackRefsStep", 'stack_refs', False),