From 05d3c960fca68340232c6cd59c67662bdae7770b Mon Sep 17 00:00:00 2001 From: Melanie Clarke Date: Tue, 17 Sep 2024 09:46:13 -0400 Subject: [PATCH] JP-3639: Implement 1/f noise correction for ramp data (#8669) --- CHANGES.rst | 16 + docs/jwst/clean_flicker_noise/arguments.rst | 71 + docs/jwst/clean_flicker_noise/index.rst | 13 + docs/jwst/clean_flicker_noise/main.rst | 253 +++ .../clean_flicker_noise/reference_files.rst | 3 + docs/jwst/nsclean/arguments.rst | 62 +- docs/jwst/nsclean/main.rst | 168 +- docs/jwst/package_index.rst | 1 + docs/jwst/pipeline/calwebb_detector1.rst | 88 +- jwst/clean_flicker_noise/__init__.py | 3 + .../clean_flicker_noise.py | 1522 +++++++++++++++++ .../clean_flicker_noise_step.py | 146 ++ jwst/{nsclean => clean_flicker_noise}/lib.py | 3 +- jwst/clean_flicker_noise/tests/__init__.py | 0 .../tests/test_clean_flicker_noise.py | 844 +++++++++ .../tests/test_clean_flicker_noise_step.py | 82 + jwst/lib/suffix.py | 2 + jwst/msaflagopen/msaflag_open.py | 21 +- jwst/msaflagopen/tests/test_msa_open.py | 96 +- jwst/nsclean/nsclean.py | 508 ------ jwst/nsclean/nsclean_step.py | 103 +- jwst/pipeline/calwebb_detector1.py | 5 + jwst/pipeline/calwebb_spec2.py | 4 +- jwst/pipeline/clean_noise.cfg | 2 + jwst/pipeline/tests/test_calwebb_spec2.py | 2 +- jwst/regtest/test_miri_image.py | 43 + jwst/regtest/test_nircam_image.py | 41 + jwst/regtest/test_niriss_image.py | 38 +- jwst/regtest/test_nirspec_irs2_detector1.py | 41 + jwst/step.py | 2 + jwst/stpipe/integration.py | 1 + 31 files changed, 3397 insertions(+), 787 deletions(-) create mode 100644 docs/jwst/clean_flicker_noise/arguments.rst create mode 100644 docs/jwst/clean_flicker_noise/index.rst create mode 100644 docs/jwst/clean_flicker_noise/main.rst create mode 100644 docs/jwst/clean_flicker_noise/reference_files.rst create mode 100644 jwst/clean_flicker_noise/__init__.py create mode 100644 jwst/clean_flicker_noise/clean_flicker_noise.py create mode 100644 jwst/clean_flicker_noise/clean_flicker_noise_step.py rename jwst/{nsclean => clean_flicker_noise}/lib.py (99%) create mode 100644 jwst/clean_flicker_noise/tests/__init__.py create mode 100644 jwst/clean_flicker_noise/tests/test_clean_flicker_noise.py create mode 100644 jwst/clean_flicker_noise/tests/test_clean_flicker_noise_step.py delete mode 100644 jwst/nsclean/nsclean.py create mode 100644 jwst/pipeline/clean_noise.cfg diff --git a/CHANGES.rst b/CHANGES.rst index b9f075c9c6..51fdfe145b 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -42,6 +42,18 @@ badpix_selfcal - Subtract pedestal dark when constructing min array across selfcal exposures for MIRI MRS data. [#8771] +calwebb_detector1 +----------------- + +- Added the optional ``clean_flicker_noise`` step between ``jump`` and + ``ramp_fit``. [#8669] + +clean_flicker_noise +------------------- + +- Implemented this new optional step to clean transient flicker noise (e.g. 1/f noise) + from group images in ramp data. [#8669] + cube_build ---------- @@ -115,6 +127,10 @@ nsclean Due to the computational costs of matrix operations, this is a large speedup that has little effect on the results. [#8745] +- Merged implementation with the new ``clean_flicker_noise`` step. This step + can still be called from the ``calwebb_spec2`` pipeline on NIRSpec rate + data, but it is now deprecated. [#8669] + outlier_detection ----------------- diff --git a/docs/jwst/clean_flicker_noise/arguments.rst b/docs/jwst/clean_flicker_noise/arguments.rst new file mode 100644 index 0000000000..8fc131709e --- /dev/null +++ b/docs/jwst/clean_flicker_noise/arguments.rst @@ -0,0 +1,71 @@ +.. _clean_flicker_noise_arguments: + +Step Arguments +============== + +The ``clean_flicker_noise`` step has the following optional arguments to control +the behavior of the processing. + +``--fit_method`` (str, default='median') + The noise fitting algorithm to use. Options are 'fft' and 'median'. + +``--fit_by_channel`` (boolean, default=False) + If set, flicker noise is fit independently for each detector channel. + Ignored for MIRI, for subarray data, and for `fit_method` = 'fft'. + +``--background_method`` (str, default='median') + If 'median', the preliminary background to remove and restore + after fitting the noise is a simple median of the background data. + If 'model', the background data is fit with a low-resolution model + via `~photutils.background.Background2D`. + If None, the background value is set to 0.0. + +``--background_box_size`` (list of int, default=(32,32)) + Box size for the data grid used by `~photutils.background.Background2D` + when `background_method` = 'model'. For best results, use a + box size that evenly divides the input image shape. + +``--mask_science_regions`` (boolean, default=False) + For NIRSpec, mask regions of the image defined by WCS bounding + boxes for slits/slices, as well as any regions known to be + affected by failed-open MSA shutters. For MIRI imaging, mask + regions of the detector not used for science. + +``--n_sigma`` (float, default=2.0) + The sigma-clipping threshold to use when searching for outliers + and illuminated pixels to be excluded from use in the background + and noise fitting processes. + +``--fit_histogram`` (boolean, default=False) + If set, the 'sigma' used with `n_sigma` for clipping outliers + is derived from a Gaussian fit to a histogram of values. + Otherwise, a simple iterative sigma clipping is performed. + +``--single_mask`` (boolean, default=True) + If set, a single mask will be created, regardless of + the number of input integrations. Otherwise, the mask will + be a 3D cube, with one plane for each integration. + +``--user_mask`` (string, default=None) + Path to a user-supplied mask file. If supplied, the mask is used + directly and the process of creating a scene mask in the step is + skipped. + + The mask file must contain either a `~jwst.datamodels.ImageModel` + or a `~jwst.datamodels.CubeModel`, with image dimensions matching + the input science data. If an ImageModel is provided, the same + mask will be used for all integrations. If a CubeModel is provided, + the number of slices must equal the number of integrations in + the input science data. + +``--save_mask`` (boolean, default=False) + If set, the mask constructed by the step will be saved to a file + with suffix 'mask'. + +``--save_background`` (boolean, default=False) + If set, the background fit to the group diff images will be saved + to a file with suffix 'flicker_bkg'. + +``--save_noise`` (boolean, default=False) + If set, the residual noise fit and removed from the input data + will be saved to a file with suffix 'flicker_noise'. diff --git a/docs/jwst/clean_flicker_noise/index.rst b/docs/jwst/clean_flicker_noise/index.rst new file mode 100644 index 0000000000..aa8b5e426d --- /dev/null +++ b/docs/jwst/clean_flicker_noise/index.rst @@ -0,0 +1,13 @@ +.. _clean_flicker_noise_step: + +=================== +Clean Flicker Noise +=================== +.. toctree:: + :maxdepth: 2 + + main.rst + reference_files.rst + arguments.rst + +.. automodapi:: jwst.clean_flicker_noise diff --git a/docs/jwst/clean_flicker_noise/main.rst b/docs/jwst/clean_flicker_noise/main.rst new file mode 100644 index 0000000000..c7356c690a --- /dev/null +++ b/docs/jwst/clean_flicker_noise/main.rst @@ -0,0 +1,253 @@ +Description +=========== + +:Class: `jwst.clean_flicker_noise.CleanFlickerNoiseStep` +:Alias: clean_flicker_noise + +Overview +-------- +The ``clean_flicker_noise`` step removes flicker noise from calibrated +ramp images, after the :ref:`jump ` step and prior to +performing the :ref:`ramp_fitting ` step. + +For NIR detectors, the noise addressed by this step is 1/f noise, which +appears as faint banding along the detector slow axis. For NIRSpec and +NIRISS, 1/f noise looks like vertical stripes; for NIRCam, it appears +as horizontal stripes. + +MIRI images also often show a faint vertical banding, similar to 1/f noise +but from a different physical source. This type of flicker noise can be +corrected by similar methods, so it is also addressed by this step. + +To correct for flicker noise, the algorithm requires that the noise +generated between one group readout and the next be isolated as much +as possible from the astronomical signal. The residual noise may then +be fit in frequency space, via an FFT, or may be characterized by a +median along the rows or columns, as appropriate for the detector. +The modeled noise is then directly subtracted from each group readout. + +The correction is available for any type of NIRSpec, NIRISS, or NIRCam +exposure. For MIRI, it is available only for imaging exposures. + +Creation of a scene mask +------------------------ +One of the key components of the correction is knowing which pixels can +be used to fit the background noise. The step builds a scene mask +on the fly from a draft rate image, generated from the input ramp data, +which is used to mark usable and unusable pixels. The mask is a 2D +Boolean array, having the same size as the image, with +pixels set to True interpreted as being OK to use. + +The process of building the mask varies somewhat depending on the +observing mode of the image being processed. Some features are common +to all modes, while others are mode-specific. The following sections +describe each type of masking that can be applied. At the end, there +is a summary of the types applied to each instrument mode. + +The user-settable step parameter `save_mask` can be used to save the +scene mask to a file, if desired (see the +:ref:`step arguments `). + +Note that a user may also supply their own mask image as input to the step, +in which case the process of creating a mask is skipped. The step parameter +`user_mask` is used to specify an input mask. If specified, the input +mask must contain a datamodel matching the shape of a 'rate' or 'rateints' +file generated for the input science data (either `~jwst.datamodels.ImageModel` +or `~jwst.datamodels.CubeModel`). To generate a custom mask, it may be +easiest to save the mask output from a default run of the ``clean_flicker_noise`` +step on the input data, then modify it as needed. + +NIRSpec IFU Slices +^^^^^^^^^^^^^^^^^^ +For IFU images, the majority of the mask is based on knowing which +pixels are contained within the footprints of the IFU slices. To do +this, the image's World Coordinate System (WCS) object is queried in +order to determine which pixels are contained within each of the 30 +slices. Pixels within each slice are set to False (do not use) in the +mask. + +NIRSpec MOS/FS Slits +^^^^^^^^^^^^^^^^^^^^ +The footprints of each open MOS slitlet or fixed slit are flagged in +a similar way as IFU slices. For MOS and FS images, the WCS object is +queried to determine which pixels are contained within each open +slit/slitlet and they are set to False in the mask. + +NIRSpec MSA Failed Open Shutters +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Pixels affected by stuck open MSA shutters are masked, because they +may contain signal. This is accomplished by setting all pixels flagged by the +:ref:`msaflagopen ` step with DQ value "MSA_FAILED_OPEN" +to False in the mask. + +NIRSpec Fixed-Slit Region Pixels +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Full-frame MOS and IFU images may contain signal from the always open +fixed slits, which appear in a fixed region in the middle of each image. +The entire region containing the fixed slits is masked out when +processing MOS and IFU images. The masked region is currently hardwired +in the step to image indexes [1:2048, 923:1116], where the indexes are +in x, y order and in 1-indexed values. + +Note, however, that it is possible to plan one or more fixed slit targets +alongside MSA slitlets in MOS observations. In this situation, the fixed +slit region is not automatically masked. + +MIRI Imaging Metering Structure +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +The MIRI imager has a metering structure covering a large part of the +detector area. These regions must be masked in order to fit and +remove the background data in the science areas of the detector. +Regions marked with DQ value "DO_NOT_USE" by the +:ref:`flat_field ` step are set to False in the +scene mask. + +Missing Data +^^^^^^^^^^^^ +Any pixel in the draft rate image that has a value of NaN or exactly zero +is flagged as False in the mask. This typically includes any reference +pixels that are present in the exposure. + +Outliers +^^^^^^^^ +For imaging modes, bright, compact sources must be distinguished +from the background and masked in order to fit and remove the +smooth background level. + +For spectral modes that already have significant masking applied, +pixels in the unilluminated areas of the region can still contain anomalous +signal, due to uncaught cosmic rays, hot pixels, etc. + +For both modes, a sigma-clipping routine is employed to find such outliers +within the draft rate image and set them to False in the mask. All pixels with +values greater than :math:`median+n\_sigma*sigma` are assumed to contain +signal and are set to False in the scene mask. In addition, all pixels +with values less than :math:`median-3.0*sigma` are assumed to be bad pixels, +and are also set to False in the scene mask. + +Mode-Specific Masking Steps +^^^^^^^^^^^^^^^^^^^^^^^^^^^ +The following table indicates which flavors of masking are applied to +images from each instrument and observing mode. + +.. |c| unicode:: U+2713 .. checkmark + ++--------------------------+-----+-----+-----+-------+--------+--------+ +| | NIRSpec | MIRI | NIRCam | NIRISS | ++--------------------------+-----+-----+-----+-------+--------+--------+ +| | IFU | MOS | FS | Image | All | All | ++==========================+=====+=====+=====+=======+========+========+ +| IFU Slices\ :sup:`1` | |c| | | | | | | ++--------------------------+-----+-----+-----+-------+--------+--------+ +| Slits/Slitlets\ :sup:`1` | | |c| | |c| | | | | ++--------------------------+-----+-----+-----+-------+--------+--------+ +| MSA_FAILED_OPEN\ :sup:`1`| |c| | |c| | | | | | ++--------------------------+-----+-----+-----+-------+--------+--------+ +| Non-science\ :sup:`1` | | | | |c| | | | ++--------------------------+-----+-----+-----+-------+--------+--------+ +| FS Region\ :sup:`1` | |c| | |c| | | | | | ++--------------------------+-----+-----+-----+-------+--------+--------+ +| Missing Data | |c| | |c| | |c| | |c| | |c| | |c| | ++--------------------------+-----+-----+-----+-------+--------+--------+ +| Outliers | |c| | |c| | |c| | |c| | |c| | |c| | ++--------------------------+-----+-----+-----+-------+--------+--------+ + +:sup:`1`\ These steps are only applied if the +:ref:`step parameter ` +`mask_science_regions` is set to True. + +Correction Algorithm +-------------------- + +The detailed process for fitting and removing flicker noise is as follows. +See the :ref:`step arguments ` for more +information on all referenced parameters. + +#. From the calibrated ramp input, make a draft rate (`single_mask` = True) + or rateints (`single_mask` = False) file. + +#. Create a scene mask from the rate data. + + #. If `mask_science_regions` is set and the input is NIRSpec data, + run :ref:`assign_wcs ` and + :ref:`msaflagopen ` on the draft rate data, + then mask any known science areas or failed-open MSA shutters. + + This will mask out regions that are likely to contain significant + astronomical signal. + + #. If `mask_science_regions` is set and the input is MIRI imaging data, + run :ref:`flat_field ` on the draft rate data, + and extract just the DQ plane from the output. Pixels flagged + as 'DO_NOT_USE' by the flat fielding process are masked. + + This will mask out regions of the detector under the metering + structure. + + #. Iteratively sigma clip the data to get a center value (mean or median) + and sigma value (standard deviation). + + #. If `fit_histogram` is set, compute a histogram from 4-sigma clipped + values and fit a Gaussian to it to refine the center and sigma values. + + #. Mask data more than 3 * sigma below the center as bad values. + + #. Mask data more than `n_sigma` * sigma above the center as signal + (not background). + +#. Iterate over each integration and group in the data, to fit and correct + for noise. + + #. Make a diff image (current group – previous group) to correct. + + #. Fit and remove a background level, using the scene mask to identify + background pixels. + + #. Clip the background data in the diff image to remove more outliers. + + #. If `background_method` = 'median', the background value is a simple + median of the remaining values. + + #. If `background_method` = 'model', the background data is fit with + a low-resolution model via the photutils + `Background2D `_ + utility. The resolution box size is set by `background_box_size`. + + #. Subtract the background level from the diff image and clip again + to `n_sigma` * sigma, with sigma recomputed from the + background-subtracted data in the remaining background pixels. + + #. Fit and remove the residual noise in the background-subtracted image. + + #. If `fit_method` = 'fft', the ``nsclean`` library is called to fit + and remove the noise in frequency space. + + #. If `fit_method` = 'median', the noise is fit with a simple median + along the appropriate detector axis and subtracted from the + background-subtracted image. + + If `fit_by_channel` = True, and the data is a NIR full-frame exposure, + the median value is computed and subtracted independently for each + detector channel. + + #. Restore the background level to the cleaned, background-subtracted + diff image. + + #. Add the cleaned diff back to a cleaned version of the previous + group image. + +References +========== + +The FFT cleaning algorithm implementation is based on NSClean, +developed by Bernard Rauscher. Details on the source of the correlated +noise and the algorithm used by the ``nsclean`` library to fit and +remove it can be found in +`Rauscher 2024 `__. + +The background fitting and median cleaning algorithm are based on +the image1overf algorithm, developed by Chris Willott, and available +on GitHub at `chriswillott/jwst `__. +The algorithm was adapted to the `clean_flicker_noise` step and is released +under the BSD license for the JWST calibration pipeline by permission +of the author. diff --git a/docs/jwst/clean_flicker_noise/reference_files.rst b/docs/jwst/clean_flicker_noise/reference_files.rst new file mode 100644 index 0000000000..552e1717a5 --- /dev/null +++ b/docs/jwst/clean_flicker_noise/reference_files.rst @@ -0,0 +1,3 @@ +Reference Files +=============== +The ``clean_flicker_noise`` step does not use any reference files. diff --git a/docs/jwst/nsclean/arguments.rst b/docs/jwst/nsclean/arguments.rst index c91bfaba3e..9772ce6497 100644 --- a/docs/jwst/nsclean/arguments.rst +++ b/docs/jwst/nsclean/arguments.rst @@ -6,19 +6,65 @@ Step Arguments The ``nsclean`` step has the following optional arguments to control the behavior of the processing. +``--fit_method`` (str, default='fft') + The noise fitting algorithm to use. Options are 'fft' and 'median'. + +``--fit_by_channel`` (boolean, default=False) + If set, flicker noise is fit independently for each detector channel. + Ignored for subarray data and for `fit_method` = 'fft'. + +``--background_method`` (str, default=None) + If 'median', the preliminary background to remove and restore + after fitting the noise is a simple median of the background data. + If 'model', the background data is fit with a low-resolution model + via `~photutils.background.Background2D`. + If None, the background value is set to 0.0. + +``--background_box_size`` (list of int, default=(32,32)) + Box size for the data grid used by `~photutils.background.Background2D` + when `background_method` = 'model'. For best results, use a + box size that evenly divides the input image shape. + ``--mask_spectral_regions`` (boolean, default=True) - Mask regions in IFU and MOS images that are within the bounding boxes - for each slice or slitlet defined in the WCS object of the image. + Mask regions of the image defined by WCS bounding + boxes for slits/slices, as well as any regions known to be + affected by failed-open MSA shutters. ``--n_sigma`` (float, default=5.0) The sigma-clipping threshold to use when searching for outliers - and illuminated pixels to be excluded from use in the fitting - process. + and illuminated pixels to be excluded from use in the background + and noise fitting processes. -``--save_mask`` (boolean, default=False) - A flag to indicate whether the mask constructed by the step - should be saved to a file. +``--fit_histogram`` (boolean, default=False) + If set, the 'sigma' used with `n_sigma` for clipping outliers + is derived from a Gaussian fit to a histogram of values. + Otherwise, a simple iterative sigma clipping is performed. + +``--single_mask`` (boolean, default=False) + If set, a single mask will be created, regardless of + the number of input integrations. Otherwise, the mask will + be a 3D cube, with one plane for each integration. ``--user_mask`` (string, default=None) Path to a user-supplied mask file. If supplied, the mask is used - directly and the process of creating a mask in the step is skipped. + directly and the process of creating a scene mask in the step is + skipped. + + The mask file must contain either a `~jwst.datamodels.ImageModel` + or a `~jwst.datamodels.CubeModel`, with image dimensions matching + the input science data. If an ImageModel is provided, the same + mask will be used for all integrations. If a CubeModel is provided, + the number of slices must equal the number of integrations in + the input science data. + +``--save_mask`` (boolean, default=False) + If set, the mask constructed by the step will be saved to a file + with suffix 'mask'. + +``--save_background`` (boolean, default=False) + If set, the background fit to the group diff images will be saved + to a file with suffix 'flicker_bkg'. + +``--save_noise`` (boolean, default=False) + If set, the residual noise fit and removed from the input data + will be saved to a file with suffix 'flicker_noise'. diff --git a/docs/jwst/nsclean/main.rst b/docs/jwst/nsclean/main.rst index 161e84d9d1..40dcda7a69 100644 --- a/docs/jwst/nsclean/main.rst +++ b/docs/jwst/nsclean/main.rst @@ -4,152 +4,26 @@ Description :Class: `jwst.nsclean.NSCleanStep` :Alias: nsclean -Overview -======== -The ``nsclean`` step applies an algorithm for removing correlated read -noise from NIRSpec images. The noise often appears as faint vertical -banding and so-called "picture frame noise." The algorithm uses dark -(unilluminated) areas of an image to fit a background model in Fourier -space. When the fit is subtracted, it removes nearly all correlated noise. -Compared to simpler strategies, like subtracting a rolling median, this -algorithm is more thorough and uniform. It is also computationally -undemanding, typically requiring only a few seconds to clean a full-frame -image. - -The correction can be applied to any type of NIRSpec exposure, including -IFU, MOS, fixed slit, and Bright Object Time Series (BOTS), in both full-frame -and subarray readouts. Time series (3D) data are corrected one integration -at a time. - -Details on the source of the correlated noise and the algorithm used -in the ``nsclean`` step to fit and remove it can be found in -`Rauscher 2023 `_. - -Upon completion of the step, the step status keyword "S_NSCLEN" gets set -to "COMPLETE" in the output science data. - -Assumptions -=========== -As described below, the creation of a pixel mask depends on the presence -of a World Coordinate System (WCS) object for the image, which is -constructed by the :ref:`assign_wcs ` step. -In addition, creating a mask for IFU and MOS images depends on -the presence of DQ flags assigned by the -:ref:`msaflagopen ` step. -It is therefore required that those steps be run before attempting to -apply ``nsclean``. - -Creation of an image mask -========================= -One of the key components of the correction is knowing which pixels are -unilluminated and hence can be used in fitting the background noise. -The step builds a mask on the fly for each image, which is used to mark -useable and unuseable pixels. The mask is a 2D boolean array, having the same -size as the image, with pixels set to True interpreted as being OK to use. -The process of building the mask varies somewhat depending on the -observing mode of the image being processed. Some features are common -to all modes, while others are mode-specific. The following sections -describe each type of masking that can be applied and at the end there -is a summary of the types applied to each image mode. - -The user-settable step parameter `save_mask` can be used to save the -mask to a file, if desired (see :ref:`nsclean step arguments `). - -Note that a user may supply their own mask image as input to the step, -in which case the process of creating a mask is skipped. The step parameter -`user_mask` is used to specify an input mask. - -IFU Slices ----------- -For IFU images the majority of the mask is based on knowing which -pixels are contained within the footprints of the IFU slices. To do -this, the image's World Coordinate System (WCS) object is queried in -order to determine which pixels are contained within each of the 30 -slices. Pixels within each slice are set to False (do not use) in the -mask. +.. note:: -MOS/FS Slits ------------- -The footprints of each open MOS slitlet or fixed slit are flagged in -a similar way as IFU slices. For MOS and FS images, the WCS object is -queried to determine which pixels are contained within each open -slit/slitlet and they are set to False in the mask. + This step, which is intended to be called in the + :ref:`calwebb_spec2 ` pipeline on NIRSpec rate data, + is now deprecated. In future builds, it will be replaced by + the :ref:`clean_flicker_noise ` + step, called in the :ref:`calwebb_detector1 ` + pipeline on ramp data. -MSA Failed Open Shutters ------------------------- -Pixels affected by stuck open MSA shutters are masked, because they -may contain signal. This is accomplished by setting all pixels flagged by the -:ref:`msaflagopen ` step with DQ value "MSA_FAILED_OPEN" -to False in the mask. - -NaN Pixels ----------- -Any pixel in the input image that has a value of NaN is temporarily reset -to zero for input to the fitting routine and flagged as False in the mask. -Upon completion of the noise subtraction, this population of pixels is -set back to NaN again in the output (corrected) image. - -Fixed-Slit Region Pixels ------------------------- -Full-frame MOS and IFU images may contain signal from the always open -fixed slits, which appear in fixed region in the middle of each image. -The entire region containing the fixed slits is masked out when -processing MOS and IFU images. The masked region is currently hardwired -in the step to image indexes [1:2048, 923:1116], where the indexes are -in x, y order and in 1-indexed values. - -Note, however, that it is possible to plan one or more fixed slit targets -alongside MSA slitlets in MOS observations. In this situation, the fixed -slit region is not automatically masked. - -Left/Right Reference Pixel Columns ----------------------------------- -Full-frame images contain 4 columns of reference pixels on the left and -right edges of the image. These are not to be used in the fitting -algorithm and hence are set to False in the mask. - -Outliers --------- -Pixels in the unilluminated areas of the region can contain anomalous -signal, due to uncaught Cosmic Rays, hot pixels, etc. A sigma-clipping -routine is employed to find such outliers within the input image and set -them to False in the mask. All pixels with values greater than -:math:`median+n_sigma*sigma` are set to False in the mask. -Here `median` and `sigma` are computed -from the image using the astropy.stats `sigma_clipped_stats` routine, -using the image mask to exclude pixels that have already been flagged -and a clipping level of 5 sigma. `n_sigma` is a user-settable step -parameter, with a default value of 5.0 -(see :ref:`nsclean step arguments `). - -Mode-Specific Masking Steps ---------------------------- -The following table indicates which flavors of masking are applied to -images from each type of observing mode. - -.. |c| unicode:: U+2713 .. checkmark - -+--------------------------+-----+-----+-----+ -| | | Mode| | -+--------------------------+-----+-----+-----+ -| Masking | IFU | MOS | FS | -+==========================+=====+=====+=====+ -| IFU Slices\ :sup:`1` | |c| | | | -+--------------------------+-----+-----+-----+ -| Slits/Slitlets\ :sup:`1` | | |c| | |c| | -+--------------------------+-----+-----+-----+ -| MSA_FAILED_OPEN | |c| | |c| | |c| | -+--------------------------+-----+-----+-----+ -| NaN Pixels | |c| | |c| | |c| | -+--------------------------+-----+-----+-----+ -| FS Region | |c| | |c| | | -+--------------------------+-----+-----+-----+ -| Reference Pix | |c| | |c| | |c| | -+--------------------------+-----+-----+-----+ -| Outliers | |c| | |c| | |c| | -+--------------------------+-----+-----+-----+ - -:sup:`1`\ The application of these steps can be turned on and off via -the step parameter `mask_spectral_regions`. This parameter controls -whether the "IFU Slices" and "Slits/Slitlets" portions of the masking -are applied. +Overview +======== +This step currently runs a version of the +:ref:`clean_flicker_noise ` algorithm, +with slightly different parameters and default values, intended +to be backwards-compatible with the previous implementation of +the step developed specifically for NIRSpec data. The only +significant difference is that in this step, the cleaning is +performed directly on the rate data, rather than group images. + +See the :ref:`clean_flicker_noise ` +documentation for more information, or the +:ref:`nsclean step arguments ` for the default +values used for this step. diff --git a/docs/jwst/package_index.rst b/docs/jwst/package_index.rst index da12f11e51..4bd3b28ba1 100644 --- a/docs/jwst/package_index.rst +++ b/docs/jwst/package_index.rst @@ -16,6 +16,7 @@ Package Index badpix_selfcal/index.rst barshadow/index.rst charge_migration/index.rst + clean_flicker_noise/index.rst combine_1d/index.rst cube_build/index.rst dark_current/index.rst diff --git a/docs/jwst/pipeline/calwebb_detector1.rst b/docs/jwst/pipeline/calwebb_detector1.rst index e66d135afb..5083b0a618 100644 --- a/docs/jwst/pipeline/calwebb_detector1.rst +++ b/docs/jwst/pipeline/calwebb_detector1.rst @@ -45,52 +45,52 @@ on either the first group or frame zero pixel values. .. |check| unicode:: U+2713 .. checkmark -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| Near-IR | MIRI | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| Step | Non-TSO | TSO | Step | Non-TSO | TSO | -+======================================================================+=========+=========+=========================================+=========+=========+ -| :ref:`group_scale ` | |check| | |check| | :ref:`group_scale ` | |check| | |check| | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| :ref:`dq_init ` | |check| | |check| | :ref:`dq_init ` | |check| | |check| | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| | | | :ref:`emicorr ` | |check| | |check| | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| :ref:`saturation ` | |check| | |check| | :ref:`saturation ` | |check| | |check| | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| :ref:`ipc ` [1]_ | | | :ref:`ipc ` | | | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| :ref:`superbias ` | |check| | |check| | :ref:`firstframe ` | |check| | | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| :ref:`refpix ` | |check| | |check| | :ref:`lastframe ` | |check| | | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| | | | :ref:`reset ` | |check| | |check| | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| :ref:`linearity ` | |check| | |check| | :ref:`linearity ` | |check| | |check| | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| :ref:`persistence ` [2]_ | |check| | | :ref:`rscd ` | |check| | | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| :ref:`dark_current ` | |check| | |check| | :ref:`dark_current ` | |check| | |check| | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| | | | :ref:`refpix ` | |check| | |check| | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| :ref:`charge_migration ` [3]_ | |check| | | | | | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| :ref:`jump ` | |check| | |check| | :ref:`jump ` | |check| | |check| | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| :ref:`ramp_fitting ` | |check| | |check| | :ref:`ramp_fitting ` | |check| | |check| | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ -| :ref:`gain_scale ` | |check| | |check| | :ref:`gain_scale ` | |check| | |check| | -+----------------------------------------------------------------------+---------+---------+-----------------------------------------+---------+---------+ - - -.. [1] By default, the parameter reference `pars-detector1pipeline` - retrieved from CRDS will skip the :ref:`ipc ` step for all instruments. ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| Near-IR | MIRI | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| Step | Non-TSO | TSO | Step | Non-TSO | TSO | ++======================================================================+=========+=========+======================================================+=========+=========+ +| :ref:`group_scale ` | |check| | |check| | :ref:`group_scale ` | |check| | |check| | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| :ref:`dq_init ` | |check| | |check| | :ref:`dq_init ` | |check| | |check| | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| | | | :ref:`emicorr ` | |check| | |check| | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| :ref:`saturation ` | |check| | |check| | :ref:`saturation ` | |check| | |check| | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| :ref:`ipc ` [1]_ | | | :ref:`ipc ` | | | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| :ref:`superbias ` | |check| | |check| | :ref:`firstframe ` | |check| | | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| :ref:`refpix ` | |check| | |check| | :ref:`lastframe ` | |check| | | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| | | | :ref:`reset ` | |check| | |check| | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| :ref:`linearity ` | |check| | |check| | :ref:`linearity ` | |check| | |check| | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| :ref:`persistence ` [2]_ | |check| | | :ref:`rscd ` | |check| | | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| :ref:`dark_current ` | |check| | |check| | :ref:`dark_current ` | |check| | |check| | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| | | | :ref:`refpix ` | |check| | |check| | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| :ref:`charge_migration ` [1]_ | |check| | | | | | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| :ref:`jump ` | |check| | |check| | :ref:`jump ` | |check| | |check| | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| :ref:`clean_flicker_noise ` [1]_ | |check| | |check| | :ref:`clean_flicker_noise `| |check| | |check| | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| :ref:`ramp_fitting ` | |check| | |check| | :ref:`ramp_fitting ` | |check| | |check| | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ +| :ref:`gain_scale ` | |check| | |check| | :ref:`gain_scale ` | |check| | |check| | ++----------------------------------------------------------------------+---------+---------+------------------------------------------------------+---------+---------+ + + +.. [1] By default, this step is skipped in the ``calwebb_detector1`` pipeline + for all instruments and modes. .. [2] The :ref:`persistence ` step is currently hardwired to be skipped in the `Detector1Pipeline` module for all NIRSpec exposures. -.. [3] By default, the :ref:`charge_migration ` step is skipped in - the `Detector1Pipeline` module for all instruments. - + Arguments --------- diff --git a/jwst/clean_flicker_noise/__init__.py b/jwst/clean_flicker_noise/__init__.py new file mode 100644 index 0000000000..e9bd9a65dc --- /dev/null +++ b/jwst/clean_flicker_noise/__init__.py @@ -0,0 +1,3 @@ +from .clean_flicker_noise_step import CleanFlickerNoiseStep + +__all__ = ['CleanFlickerNoiseStep'] diff --git a/jwst/clean_flicker_noise/clean_flicker_noise.py b/jwst/clean_flicker_noise/clean_flicker_noise.py new file mode 100644 index 0000000000..2a41da373c --- /dev/null +++ b/jwst/clean_flicker_noise/clean_flicker_noise.py @@ -0,0 +1,1522 @@ +import logging +import warnings + +import gwcs +import numpy as np +from astropy.stats import sigma_clipped_stats, SigmaClip +from astropy.utils.exceptions import AstropyUserWarning +from photutils.background import Background2D, MedianBackground +from scipy.optimize import curve_fit +from stdatamodels.jwst.datamodels import dqflags + +from jwst import datamodels +from jwst.assign_wcs import nirspec, AssignWcsStep +from jwst.flatfield import FlatFieldStep +from jwst.clean_flicker_noise.lib import NSClean, NSCleanSubarray +from jwst.lib.basic_utils import LoggingContext +from jwst.msaflagopen import MSAFlagOpenStep +from jwst.ramp_fitting import RampFitStep + + +log = logging.getLogger(__name__) +log.setLevel(logging.DEBUG) + +# Fixed slit region to mask, for NIRSpec MOS and IFU data +# Values are y start and stop indices, for the edges of the +# region to mask. +NRS_FS_REGION = [922, 1116] + + +def make_rate(input_model, input_dir='', return_cube=False): + """ + Make a rate model from a ramp model. + + Parameters + ---------- + input_model : `~jwst.datamodel.RampModel` + Input ramp model. + + input_dir : str + Path to the input directory. Used by sub-steps (e.g. assign_wcs + for NIRSpec MOS data) to find auxiliary data. + + return_cube : bool, optional + If set, a CubeModel will be returned, with a separate + rate for each integration. Otherwise, an ImageModel is returned + with the combined rate for the full observation. + + Returns + ------- + rate_model : `~jwst.datamodel.ImageModel` or `~jwst.datamodel.CubeModel` + The rate or rateints model. + """ + # Call the ramp fit step on a copy of the input + # Use software default values for parameters + + log.info("Creating draft rate file for scene masking") + step = RampFitStep() + step.input_dir = input_dir + with LoggingContext(step.log, level=logging.WARNING): + # Note: the copy is currently needed because ramp fit + # closes the input model when it's done, and we need + # it to stay open. + rate, rateints = step.run(input_model.copy()) + + if return_cube: + output_model = rateints + rate.close() + del rate + else: + output_model = rate + rateints.close() + del rateints + + return output_model + + +def post_process_rate(input_model, input_dir='', assign_wcs=False, + msaflagopen=False, flat_dq=False): + """ + Post-process the input rate model, as needed. + + Parameters + ---------- + input_model : `~jwst.datamodel.ImageModel` or `~jwst.datamodel.CubeModel` + Input rate model. + + input_dir : str + Path to the input directory. Used by sub-steps (e.g. assign_wcs + for NIRSpec MOS data) to find auxiliary data. + + assign_wcs : bool, optional + If set and the input does not already have a WCS assigned, + the assign_wcs step will be called on the rate model. + + msaflagopen : bool, optional + If set, the msaflagopen step will be called on the rate model. + If a WCS is not already present, assign_wcs will be called first. + + flat_dq : bool, optional + If set, the flat_field step will be run on the input model. DQ + flags are retrieved from the output and added to the input + model's DQ array. The rate data is not modified. + + Returns + ------- + output_model : `~jwst.datamodel.ImageModel` or `~jwst.datamodel.CubeModel` + The updated model. + """ + output_model = input_model + + # If needed, assign a WCS + if (assign_wcs or msaflagopen) and not hasattr(output_model.meta, 'wcs'): + log.info("Assigning a WCS for scene masking") + step = AssignWcsStep() + step.input_dir = input_dir + with LoggingContext(step.log, level=logging.WARNING): + output_model = step.run(output_model) + + # If needed, flag open MSA shutters + if msaflagopen: + log.info("Flagging failed-open MSA shutters for scene masking") + step = MSAFlagOpenStep() + step.input_dir = input_dir + with LoggingContext(step.log, level=logging.WARNING): + output_model = step.run(output_model) + + # If needed, draft a flat correction to retrieve non-science areas + if flat_dq: + log.info("Retrieving flat DQ values for scene masking") + step = FlatFieldStep() + step.input_dir = input_dir + with LoggingContext(step.log, level=logging.WARNING): + flat_corrected_model = step.run(output_model) + + # Copy out the flat DQ plane, leave the data as is + output_model.dq |= flat_corrected_model.dq + flat_corrected_model.close() + del flat_corrected_model + + return output_model + + +def mask_ifu_slices(input_model, mask): + """ + Flag pixels within IFU slices. + + Find pixels located within IFU slices, according to the WCS, + and flag them in the mask, so that they do not get used. + + Parameters + ---------- + input_model : `~jwst.datamodel.JwstDataModel` + Science data model + + mask : array-like of bool + 2D input mask that will be updated. True indicates background + pixels to be used. IFU slice regions will be set to False. + + Returns + ------- + mask : array-like of bool + 2D output mask with additional flags for science pixels + """ + log.info("Finding slice pixels for an IFU image") + + # Initialize global DQ map to all zero (OK to use) + dqmap = np.zeros_like(input_model.dq) + + # Get the wcs objects for all IFU slices + list_of_wcs = nirspec.nrs_ifu_wcs(input_model) + + # Loop over the IFU slices, finding the valid region for each + for (k, ifu_wcs) in enumerate(list_of_wcs): + + # Construct array indexes for pixels in this slice + x, y = gwcs.wcstools.grid_from_bounding_box( + ifu_wcs.bounding_box, step=(1, 1), center=True) + + # Get the world coords for all pixels in this slice; + # all we actually need are wavelengths + coords = ifu_wcs(x, y) + dq = dqmap[..., y.astype(int), x.astype(int)] + wl = coords[2] + + # set non-NaN wavelength locations as do not use (one) + valid = ~np.isnan(wl) + dq = dq[..., valid] + x = x[..., valid] + y = y[..., valid] + dq[:] = 1 + + # Copy DQ for this slice into global DQ map + dqmap[..., y.astype(int), x.astype(int)] = dq + + # Now set all non-zero locations in the mask to False (do not use) + mask[dqmap == 1] = False + + return mask + + +def mask_slits(input_model, mask): + """ + Flag pixels within science regions. + + Find pixels located within MOS or fixed slit footprints + and flag them in the mask, so that they do not get used. + + Parameters + ---------- + input_model : `~jwst.datamodel.JwstDataModel` + Science data model. + + mask : array-like of bool + 2D input mask that will be updated. True indicates background + pixels to be used. Slit regions will be set to False. + + Returns + ------- + mask : array-like of bool + 2D output mask with additional flags for slit pixels + """ + + from jwst.extract_2d.nirspec import offset_wcs + + log.info("Finding slit/slitlet pixels") + + # Get the slit-to-msa frame transform from the WCS object + slit2msa = input_model.meta.wcs.get_transform('slit_frame', 'msa_frame') + + # Loop over the slits, marking all the pixels within each bounding + # box as False (do not use) in the mask. + # Note that for 3D masks (TSO mode), all planes will be set to the same value. + for slit in slit2msa.slits: + slit_wcs = nirspec.nrs_wcs_set_input(input_model, slit.name) + xlo, xhi, ylo, yhi = offset_wcs(slit_wcs) + mask[..., ylo:yhi, xlo:xhi] = False + + return mask + + +def clip_to_background(image, mask, sigma_lower=3.0, sigma_upper=2.0, + fit_histogram=False, lower_half_only=False, + verbose=False): + """ + Flag signal and bad pixels in the image mask. + + Given an image, estimate the background level and a sigma value for the + mean background. + + The center and sigma may be calculated from a simple sigma-clipped + median and standard deviation, or may be refined by fitting a Gaussian + to a histogram of the values. In that case, the center is the + Gaussian mean and the sigma is the Gaussian width. + + Pixels above the center + sigma_upper * sigma are assumed to + have signal; those below this level are assumed to be background pixels. + + Pixels less than center - sigma_lower * sigma are also excluded as bad values. + + The input mask is updated in place. + + Parameters + ---------- + image : array-like of float + 2D image containing signal and background values. + + mask : array-like of bool + 2D input mask to be updated. True indicates background + pixels to be used. Regions containing signal or outliers + will be set to False. + + sigma_lower : float, optional + The number of standard deviations to use as the lower bound + for the clipping limit. Values below this limit are marked + False in the mask. + + sigma_upper : float, optional + The number of standard deviations to use as the upper bound + for the clipping limit. Values above this limit are marked + False in the mask. + + fit_histogram : bool, optional + If set, the center value and standard deviation used with + `sigma_lower` and `sigma_upper` for clipping outliers is derived + from a Gaussian fit to a histogram of values. Otherwise, the + center and standard deviation are derived from a simple iterative + sigma clipping. + + lower_half_only : bool, optional + If set, the data used to compute the center and standard deviation + for clipping is the lower half of the distribution only. Values + below the median are mirrored around the median value to simulate + a symmetric distribution. This is intended to account for + asymmetrical value distributions, with long tails in the upper + half of the distribution, due to diffuse emission, for example. + + verbose : bool, optional + If set, DEBUG level messages are issued with details on the + computed statistics. + """ + # Sigma limit for basic stats + sigma_limit = 3.0 + + # Check mask for any valid data before proceeding + if not np.any(mask): + return + + # Initial iterative sigma clip + with warnings.catch_warnings(): + warnings.filterwarnings( + action="ignore", category=AstropyUserWarning) + mean, median, sigma = sigma_clipped_stats( + image, mask=~mask, sigma=sigma_limit) + if fit_histogram: + center = mean + else: + center = median + if verbose: + log.debug('From initial sigma clip:') + log.debug(f' center: {center:.5g}') + log.debug(f' sigma: {sigma:.5g}') + + # If desired, use only the lower half of the data distribution + if lower_half_only: + lower_half_idx = mask & (image <= center) + data_for_stats = np.concatenate( + ((image[lower_half_idx] - center), + (center - image[lower_half_idx]))) + center + + # Redo stats on lower half of distribution + with warnings.catch_warnings(): + warnings.filterwarnings( + action="ignore", category=AstropyUserWarning) + mean, median, sigma = sigma_clipped_stats( + data_for_stats, sigma=sigma_limit) + if fit_histogram: + center = mean + else: + center = median + if verbose: + log.debug('From lower half distribution:') + log.debug(f' center: {center:.5g}') + log.debug(f' sigma: {sigma:.5g}') + else: + data_for_stats = image[mask] + + # Refine sigma and center from a fit to a histogram, if desired + if fit_histogram: + try: + hist, edges = np.histogram( + data_for_stats, bins=2000, + range=(center - 4. * sigma, center + 4. * sigma)) + except ValueError: + log.error('Histogram failed; using clip center and sigma.') + hist, edges = None, None + + param_opt = None + if hist is not None: + values = (edges[1:] + edges[0:-1]) / 2. + ind = np.argmax(hist) + mode_estimate = values[ind] + + # Fit a Gaussian profile to the histogram + def gaussian(x, g_amp, g_mean, g_sigma): + return g_amp * np.exp(-0.5 * ((x - g_mean) / g_sigma) ** 2) + + param_start = (hist[ind], mode_estimate, sigma) + bounds = [(0, values[0], 0), + (np.inf, values[-1], values[-1] - values[0])] + try: + param_opt, _ = curve_fit(gaussian, values, hist, p0=param_start, + bounds=bounds) + except RuntimeError: + log.error('Gaussian fit failed; using clip center and sigma.') + param_opt = None + + if verbose: + log.debug('From histogram:') + log.debug(f' mode estimate: {mode_estimate:.5g}') + log.debug(f' range of values in histogram: ' + f'{values[0]:.5g} to {values[-1]:.5g}') + + if verbose: + log.debug('Gaussian fit results:') + if param_opt is None: + if verbose: + log.debug(' (fit failed)') + else: + if verbose: + log.debug(f' peak: {param_opt[0]:.5g}') + log.debug(f' center: {param_opt[1]:.5g}') + log.debug(f' sigma: {param_opt[2]:.5g}') + center = param_opt[1] + sigma = param_opt[2] + + # Set limits from center and sigma + background_lower_limit = center - sigma_lower * sigma + background_upper_limit = center + sigma_upper * sigma + if verbose: + log.debug(f'Mask limits: {background_lower_limit:.5g} ' + f'to {background_upper_limit:.5g}') + + # Clip bad values + bad_values = image < background_lower_limit + mask[bad_values] = False + + # Clip signal (> N sigma) + signal = image > background_upper_limit + mask[signal] = False + + +def create_mask(input_model, mask_science_regions=False, + n_sigma=2.0, fit_histogram=False, single_mask=False): + """ + Create a mask identifying background pixels. + + Parameters + ---------- + input_model : `~jwst.datamodel.JwstDataModel` + Science data model, containing rate data with all necessary + pre-processing already performed. + + mask_science_regions : bool, optional + For NIRSpec, mask regions of the image defined by WCS bounding + boxes for slits/slices, as well as any regions known to be + affected by failed-open MSA shutters. This requires the + `assign_wcs` and `msaflagopen` steps to have been run on the + input_model. + For MIRI imaging, mask regions of the detector not used for science. + This requires that DO_NOT_USE flags are set in the DQ array + for the input_model. + + n_sigma : float, optional + Sigma threshold for masking outliers. + + fit_histogram : bool, optional + If set, the 'sigma' used with `n_sigma` for clipping outliers + is derived from a Gaussian fit to a histogram of values. + Otherwise, a simple iterative sigma clipping is performed. + + single_mask : bool, optional + If set, a single mask will be created, regardless of + the number of input integrations. Otherwise, the mask will + be a 3D cube, with one plane for each integration. + + Returns + ------- + mask : array-like of bool + 2D or 3D image mask + """ + exptype = input_model.meta.exposure.type.lower() + + # Initialize mask to all True. Subsequent operations will mask + # out pixels that contain signal. + mask = np.full(input_model.dq.shape, True) + + # If IFU, mask all pixels contained in the IFU slices + if exptype == 'nrs_ifu' and mask_science_regions: + mask = mask_ifu_slices(input_model, mask) + + # If MOS or FS, mask all pixels affected by open slitlets + if (exptype in ['nrs_fixedslit', 'nrs_brightobj', 'nrs_msaspec'] + and mask_science_regions): + mask = mask_slits(input_model, mask) + + # If IFU or MOS, mask pixels affected by failed-open shutters + if mask_science_regions and exptype in ['nrs_ifu', 'nrs_msaspec']: + open_pix = input_model.dq & dqflags.pixel['MSA_FAILED_OPEN'] + mask[open_pix > 0] = False + + # If MIRI imaging, mask the non-science regions: + # they contain irrelevant emission + if mask_science_regions and exptype in ['mir_image']: + non_science = (input_model.dq & dqflags.pixel['DO_NOT_USE']) > 0 + mask[non_science] = False + + # Mask any NaN pixels or exactly zero value pixels + no_data_pix = np.isnan(input_model.data) | (input_model.data == 0) + mask[no_data_pix] = False + + # If IFU or MOS, mask the fixed-slit area of the image; uses hardwired indexes + if mask_science_regions: + if exptype == 'nrs_ifu': + log.info("Masking the fixed slit region for IFU data.") + mask[..., NRS_FS_REGION[0]:NRS_FS_REGION[1], :] = False + elif exptype == 'nrs_msaspec': + # check for any slits defined in the fixed slit quadrant: + # if there is nothing there of interest, mask the whole FS region + try: + slit2msa = input_model.meta.wcs.get_transform('slit_frame', 'msa_frame') + is_fs = [s.quadrant == 5 for s in slit2msa.slits] + except (AttributeError, ValueError, TypeError): + log.warning('Slit to MSA transform not found.') + is_fs = [False] + if not any(is_fs): + log.info("Masking the fixed slit region for MOS data.") + mask[..., NRS_FS_REGION[0]:NRS_FS_REGION[1], :] = False + else: + log.info("Fixed slits found in MSA definition; " + "not masking the fixed slit region for MOS data.") + + # Mask outliers and signal using sigma clipping stats. + # For 3D data, loop over each integration separately. + if input_model.data.ndim == 3: + for i in range(input_model.data.shape[0]): + clip_to_background( + input_model.data[i], mask[i], + sigma_upper=n_sigma, fit_histogram=fit_histogram, verbose=True) + else: + clip_to_background( + input_model.data, mask, + sigma_upper=n_sigma, fit_histogram=fit_histogram, verbose=True) + + # Reduce the mask to a single plane if needed + if single_mask and mask.ndim == 3: + mask = np.all(mask, axis=0) + + # Return the mask + return mask + + +def background_level(image, mask, background_method='median', + background_box_size=None): + """ + Fit a low-resolution background level. + + Parameters + ---------- + image : array-like of float + The 2D image containing the background to fit. + + mask : array-like of bool + The mask that indicates which pixels are to be used in fitting. + True indicates a background pixel. + + background_method : {'median', 'model', None}, optional + If 'median', the preliminary background to remove and restore + is a simple median of the background data. If 'model', the + background data is fit with a low-resolution model via + `~photutils.background.Background2D`. If None, the background + value is 0.0. + + background_box_size : tuple of int, optional + Box size for the data grid used by `Background2D` when + `background_method` = 'model'. For best results, use a box size + that evenly divides the input image shape. Defaults to 32x32 if + not provided. + + Returns + ------- + background : float or array-like of float + The background level: a single value, if `background_method` + is 'median' or None, or an array matching the input image size + if `background_method` is 'model'. + """ + if background_method is None: + background = 0.0 + + else: + # Sigma limit for basic stats + sigma_limit = 3.0 + + # Flag more signal in the background subtracted image, + # with sigma set by the lower half of the distribution only + clip_to_background( + image, mask, sigma_lower=sigma_limit, sigma_upper=sigma_limit, + lower_half_only=True) + + if background_method == 'model': + sigma_clip_for_bkg = SigmaClip(sigma=sigma_limit, maxiters=5) + bkg_estimator = MedianBackground() + + if background_box_size is None: + background_box_size = [32, 32] + box_division_remainder = (image.shape[0] % background_box_size[0], + image.shape[1] % background_box_size[1]) + if not np.allclose(box_division_remainder, 0): + log.warning(f'Background box size {background_box_size} ' + f'does not divide evenly into the image ' + f'shape {image.shape}.') + + try: + with warnings.catch_warnings(): + warnings.filterwarnings( + action="ignore", category=AstropyUserWarning) + bkg = Background2D( + image, box_size=background_box_size, filter_size=(5, 5), + mask=~mask, sigma_clip=sigma_clip_for_bkg, + bkg_estimator=bkg_estimator) + background = bkg.background + except ValueError: + log.error('Background fit failed, using median value.') + background = np.nanmedian(image[mask]) + else: + background = np.nanmedian(image[mask]) + return background + + +def fft_clean_full_frame(image, mask, detector): + """ + Fit and remove background noise in frequency space for a full-frame image. + + Parameters + ---------- + image : array-like of float + The image to be cleaned. + + mask : array-like of bool + The mask that indicates which pixels are to be used in fitting. + + detector : str + The name of the detector from which the data originate. + + Returns + ------- + cleaned_image : array-like of float + The cleaned image. + """ + + # Instantiate the cleaner + cleaner = NSClean(detector, mask) + + # Clean the image + try: + cleaned_image = cleaner.clean(image, buff=True) + except np.linalg.LinAlgError: + log.warning("Error cleaning image; step will be skipped") + return None + + return cleaned_image + + +def fft_clean_subarray(image, mask, detector, npix_iter=512, + fc=(1061, 1211, 49943, 49957), + exclude_outliers=True, sigrej=4, minfrac=0.05): + """ + Fit and remove background noise in frequency space for a subarray image. + + Parameters + ---------- + image : array-like of float + The image to be cleaned. + + mask : array-like of bool + The mask that indicates which pixels are to be used in fitting. + + detector : str + The name of the detector from which the data originate. + + npix_iter : int + Number of pixels to process simultaneously. Default 512. Should + be at least a few hundred to access sub-kHz frequencies in areas + where most pixels are available for fitting. Previous default + behavior corresponds to npix_iter of infinity. + + fc : tuple + Apodizing filter definition. These parameters are tunable. The + defaults happen to work well for NIRSpec BOTS exposures. + 1) Unity gain for f < fc[0] + 2) Cosine roll-off from fc[0] to fc[1] + 3) Zero gain from fc[1] to fc[2] + 4) Cosine roll-on from fc[2] to fc[3] + Default (1061, 1211, 49943, 49957) + + exclude_outliers : bool + Find and mask outliers in the fit? Default True + + sigrej : float + Number of sigma to clip when identifying outliers. Default 4. + + minfrac : float + Minimum fraction of pixels locally available in the mask in + order to attempt a correction. Default 0.05 (i.e., 5%). + + Returns + ------- + cleaned_image : array-like of float + The cleaned image. + """ + # Flip the image to detector coords. NRS1 requires a transpose + # of the axes, while NRS2 requires a transpose and flip. + if detector == "NRS1": + image = image.transpose() + mask = mask.transpose() + else: + image = image.transpose()[::-1] + mask = mask.transpose()[::-1] + + # We must do the masking of discrepant pixels here: it just + # doesn't work if we wait and do it in the cleaner. This is + # basically copied from lib.py. Use a robust estimator for + # standard deviation, then exclude discrepant pixels and their + # four nearest neighbors from the fit. + + if exclude_outliers: + med = np.median(image[mask]) + std = 1.4825 * np.median(np.abs((image - med)[mask])) + outlier = mask & (np.abs(image - med) > sigrej * std) + + mask = mask & (~outlier) + + # also get four nearest neighbors of flagged pixels + mask[1:] = mask[1:] & (~outlier[:-1]) + mask[:-1] = mask[:-1] & (~outlier[1:]) + mask[:, 1:] = mask[:, 1:] & (~outlier[:, :-1]) + mask[:, :-1] = mask[:, :-1] & (~outlier[:, 1:]) + + # Used to determine the fitting intervals along the slow scan + # direction. Pre-pend a zero so that sum_mask[i] is equal + # to np.sum(mask[:i], axis=1). + + sum_mask = np.array([0] + list(np.cumsum(np.sum(mask, axis=1)))) + + # i1 will be the first row with a nonzero element in the mask + # imax will be the last row with a nonzero element in the mask + + nonzero_mask_element = np.sum(mask, axis=1) > 0 + + if np.sum(nonzero_mask_element) == 0: + log.warning("No good pixels in mask; step will be skipped") + return None + + i1 = np.amin(np.arange(mask.shape[0])[nonzero_mask_element]) + imax = np.amax(np.arange(mask.shape[0])[nonzero_mask_element]) + + i1_vals = [] + di_list = [] + models = [] + while i1 <= imax: + + # Want npix_iter available pixels in this section. If + # there are fewer than 1.5*npix_iter available pixels in + # the rest of the image, just go to the end. + k = 0 + for k in range(i1 + 1, imax + 2): + if (sum_mask[k] - sum_mask[i1] > npix_iter + and sum_mask[-1] - sum_mask[i1] > 1.5 * npix_iter): + break + + di = k - i1 + + i1_vals += [i1] + di_list += [di] + + # Fit this section only if at least minpct% of the pixels + # are available for finding the background. Don't flag + # outliers section-by-section; we have to do that earlier + # over the full array to get reliable values for the mean + # and standard deviation. + + if np.mean(mask[i1:i1 + di]) > minfrac: + cleaner = NSCleanSubarray(image[i1:i1 + di], mask[i1:i1 + di], + fc=fc, exclude_outliers=False) + try: + models += [cleaner.clean(return_model=True)] + except np.linalg.LinAlgError: + log.warning("Error cleaning image; step will be skipped") + return None + else: + log.warning("Insufficient reference pixels for NSClean around " + "row %d; no correction will be made here." % i1) + models += [np.zeros(image[i1:i1 + di].shape)] + + # If we have reached the end of the array, we are finished. + if k == imax + 1: + break + + # Step forward by half an interval so that we have + # overlapping fitting regions. + + i1 += max(int(np.round(di/2)), 1) + + model = np.zeros(image.shape) + tot_wgt = np.zeros(image.shape) + + # When we combine different corrections computed over + # different intervals, each one the highest weight towards the + # center of its interval and less weight towards the edge. + # Use nonzero weights everywhere so that if only one + # correction is available it gets unit weight when we + # normalize. + + for i in range(len(models)): + wgt = 1.001 - np.abs(np.linspace(-1, 1, di_list[i]))[:, np.newaxis] + model[i1_vals[i]:i1_vals[i] + di_list[i]] += wgt*models[i] + tot_wgt[i1_vals[i]:i1_vals[i] + di_list[i]] += wgt + + # don't divide by zero + tot_wgt[model == 0] = 1 + model /= tot_wgt + cleaned_image = image - model + + # Restore the cleaned image to the science frame + if detector == "NRS1": + cleaned_image = cleaned_image.transpose() + else: + cleaned_image = cleaned_image[::-1].transpose() + + return cleaned_image + + +def median_clean(image, mask, axis_to_correct, fit_by_channel=False): + """ + Fit and remove background noise via median values along one image axis. + + Parameters + ---------- + image : array-like of float + The image to be cleaned. + + mask : array-like of bool + The mask that indicates which pixels are to be used in fitting. + True indicates a background pixel. + + axis_to_correct : int + For NIR detectors, the axis to correct should be the detector slow + readout direction. Values expected are 1 or 2, following the JWST + datamodel definition (meta.subarray.slowaxis). For MIRI, flicker + noise appears along the vertical direction, so `axis_to_correct` + should be set to 1 (median along the y-axis). + + fit_by_channel : bool, optional + If set, flicker noise is fit independently for each detector channel. + Ignored for MIRI and for subarray data. + + Returns + ------- + cleaned_image : array-like of float + The cleaned image. + """ + # Masked median along slow axis + masked_image = np.ma.array(image, mask=~mask) + corrected_image = image.copy() + array_axis = axis_to_correct - 1 + + # If desired, take the median over each channel separately. + + # Full frame for the NIR detectors is 2048 x 2048, and there + # are 4 channels per detector, divided along the slow axis. + # The fit_by_channel option should only be set for full frame + # NIR data, but make sure the slow axis has the right size + # here anyway. + if fit_by_channel and image.shape[array_axis] == 2048: + n_output = 4 + channel_size = 512 + else: + # For MIRI data and subarrays, the whole image should + # be considered at once, as a single "channel" + n_output = 1 + channel_size = image.shape[array_axis] + + # Compute stripes from the median over the background data + # in the channel + cstart = 0 + cstop = channel_size + for channel in range(n_output): + if array_axis == 1: + channel_image = masked_image[:, cstart:cstop] + else: + channel_image = masked_image[cstart:cstop, :] + stripes = np.ma.median(channel_image, axis=array_axis, keepdims=True) + stripes = np.ma.filled(stripes, fill_value=0.0) + + # Remove median stripes + if array_axis == 1: + corrected_image[:, cstart:cstop] = image[:, cstart:cstop] - stripes + else: + corrected_image[cstart:cstop, :] = image[cstart:cstop, :] - stripes + + cstart += channel_size + cstop += channel_size + + return corrected_image + + +def _check_input(exp_type, fit_method): + """ + Check for valid input data and options. + + Parameters + ---------- + exp_type : str + Exposure type for the input. + fit_method : str + Noise fitting method. + + Returns + ------- + valid : bool + True if the input is valid. + """ + message = None + miri_allowed = ['MIR_IMAGE'] + if exp_type.startswith('MIR') and exp_type not in miri_allowed: + message = f"EXP_TYPE {exp_type} is not supported" + + nsclean_allowed = ['NRS_MSASPEC', 'NRS_IFU', 'NRS_FIXEDSLIT', 'NRS_BRIGHTOBJ'] + if fit_method == 'fft': + if exp_type not in nsclean_allowed: + message = f"Fit method 'fft' cannot be applied to exp_type {exp_type}" + + if message is not None: + log.warning(message) + log.warning("Step will be skipped") + return False + + return True + + +def _make_intermediate_model(input_model, intermediate_data): + """ + Make a data model to contain intermediate outputs. + + The output model type depends on the shape of the input + intermediate data. + + Parameters + ---------- + input_model : `~jwst.datamodel.JwstDataModel` + The input data. + intermediate_data : array-like + The intermediate data to save. + + Returns + ------- + intermediate_model : ~jwst.datamodel.JwstDataModel` + A model containing only the intermediate data and top-level + metadata matching the input. + """ + if intermediate_data.ndim == 4: + intermediate_model = datamodels.RampModel(data=intermediate_data) + elif intermediate_data.ndim == 3: + intermediate_model = datamodels.CubeModel(data=intermediate_data) + else: + intermediate_model = datamodels.ImageModel(data=intermediate_data) + + # Copy metadata from input model + intermediate_model.update(input_model) + return intermediate_model + + +def _standardize_parameters( + exp_type, subarray, slowaxis, background_method, fit_by_channel): + """ + Standardize input parameters. + + Check input parameters against input exposure type and assemble + values needed for subsequent correction. + + Parameters + ---------- + exp_type : str + Exposure type for the input data. + subarray : str + Subarray name for the input data. + slowaxis : int + Detector slow axis. + background_method : str + Input option for background method. + fit_by_channel : bool + Input option to fit noise by channel. + + Returns + ------- + axis_to_correct : int + Axis along which flicker noise appears for the input + exposure type. + background_method : str or None + Standardized parameter value for background correction. + fit_by_channel : bool + Standardized parameter value for the input subarray and + exposure type. + fc : tuple of int + Frequency cutoff values for use with FFT correction, + by input subarray. + """ + # Get axis to correct, by instrument + if exp_type.startswith('MIR'): + # MIRI doesn't have 1/f-noise, but it does have a vertical flickering. + # Set the axis for median correction to the y-axis. + axis_to_correct = 1 + else: + # For NIR detectors, the axis for median correction is the + # detector slowaxis. + axis_to_correct = abs(slowaxis) + + # Standardize background arguments + if str(background_method).lower() == 'none': + background_method = None + + # Check for fit_by_channel argument, and use only if data is full frame + if fit_by_channel and (subarray != 'FULL' or exp_type.startswith('MIR')): + log.warning("Fit by channel can only be used for full-frame NIR data.") + log.warning("Setting fit_by_channel to False.") + fit_by_channel = False + + # For a fft correction, ALLSLITS exposures need different + # ranges of 1/f frequencies. Be less aggressive with + # fitting higher frequencies in ALLSLITS mode. + if subarray == "ALLSLITS": + fc = (150, 200, 49943, 49957) + else: + fc = (1061, 1211, 49943, 49957) + + return axis_to_correct, background_method, fit_by_channel, fc + + +def _make_processed_rate_image( + input_model, single_mask, input_dir, exp_type, mask_science_regions): + """ + Make a draft rate image and postprocess if needed. + + Parameters + ---------- + input_model : `~jwst.datamodel.JwstDataModel` + The input data. + single_mask : bool + If set, a single scene mask is desired, so create + a single rate image. Otherwise, create a rateints + cube, with independent rate information for each + integration. + input_dir : str + Path to the input directory. Used by sub-steps to find auxiliary + data. + exp_type : str + Exposure type for the input data, used to determine + which kinds of postprocessing is necessary. + mask_science_regions : bool + If set, science regions should be masked, so run `assign_wcs` + and `msaflagopen` for NIRSpec and `flatfield` for MIRI + after creating the draft rate image. + + Returns + ------- + image_model : `~jwst.datamodel.JwstDataModel` + The processed rate image or cube. + """ + if isinstance(input_model, datamodels.RampModel): + image_model = make_rate(input_model, return_cube=(not single_mask), + input_dir=input_dir) + else: + # input is already a rate file + image_model = input_model + + # If needed, assign a WCS to the rate file, + # flag open MSA shutters, or retrieve flat DQ flags + assign_wcs = exp_type.startswith('NRS') + flag_open = (exp_type in ['NRS_IFU', 'NRS_MSASPEC']) + flat_dq = (exp_type in ['MIR_IMAGE']) + if mask_science_regions: + image_model = post_process_rate( + image_model, assign_wcs=assign_wcs, + msaflagopen=flag_open, flat_dq=flat_dq, + input_dir=input_dir) + + return image_model + + +def _make_scene_mask( + user_mask, image_model, mask_science_regions, + n_sigma, fit_histogram, single_mask, save_mask): + """ + Make a scene mask from user input or rate image. + + If provided, the user mask is opened as a datamodel and directly + returned. Otherwise, the mask is generated from the rate data in + `image_model`. + + Parameters + ---------- + user_mask : str or None + Path to user-supplied mask image. + image_model : `~jwst.datamodel.JwstDataModel` + A rate image or cube, processed as needed. + mask_science_regions: bool + For NIRSpec, mask regions of the image defined by WCS bounding + boxes for slits/slices, as well as any regions known to be + affected by failed-open MSA shutters. For MIRI imaging, mask + regions of the detector not used for science. + n_sigma : float + N-sigma rejection level for finding outliers. + fit_histogram : bool + If set, the 'sigma' used with `n_sigma` for clipping outliers + is derived from a Gaussian fit to a histogram of values. + Otherwise, a simple iterative sigma clipping is performed. + single_mask : bool + If set, a single mask will be created, regardless of + the number of input integrations. Otherwise, the mask will + be a 3D cube, with one plane for each integration. + save_mask : bool + If set, a mask model is created and returned along with the mask + array. If not, the `mask_model` returned is None. + + Returns + ------- + background_mask : array-like of bool + Mask array, with True indicating background pixels, False + indicating source pixels. + mask_model : `~jwst.datamodel.JwstDataModel` or None + A datamodel containing the background mask, if `save_mask` + is True. + """ + # Check for a user-supplied mask image. If provided, use it. + if user_mask is not None: + mask_model = datamodels.open(user_mask) + background_mask = (mask_model.data.copy()).astype(np.bool_) + mask_model.close() + del mask_model + else: + # Create the pixel mask that will be used to indicate which pixels + # to include in the 1/f noise measurements. Basically, we're setting + # all illuminated pixels to False, so that they do not get used, and + # setting all unilluminated pixels to True (so they DO get used). + log.info("Creating mask") + background_mask = create_mask( + image_model, + mask_science_regions=mask_science_regions, + n_sigma=n_sigma, fit_histogram=fit_histogram, + single_mask=single_mask) + + # Store the mask image in a model, if requested + if save_mask: + mask_model = _make_intermediate_model(image_model, background_mask) + else: + mask_model = None + + return background_mask, mask_model + + +def _check_data_shapes(input_model, background_mask): + """ + Check data shape for input model and background mask. + + Parameters + ---------- + input_model : `~jwst.datamodel.JwstDataModel` + The input data model. + background_mask : array-like of bool + The background mask. + + Returns + ------- + mismatch : bool + If True, the background mask does not match the data + and the step should be skipped. + ndim : int + Number of dimensions in the input data. + nints : int + Number of integrations in the input data. + ngroups : int + Number of groups in the input data. + """ + mismatch = True + image_shape = input_model.data.shape[-2:] + ndim = input_model.data.ndim + if ndim == 2: + nints = 1 + ngroups = 1 + else: + nints = input_model.data.shape[0] + if ndim == 3: + ngroups = 1 + else: + ngroups = input_model.data.shape[1] - 1 + + # Check for 3D mask + if background_mask.ndim == 2: + if nints > 1: + log.info("Data has multiple integrations, but mask is 2D: " + "the same mask will be used for all integrations.") + elif background_mask.shape[0] != nints: + log.warning("Mask does not match data shape. Step will be skipped.") + return mismatch, ndim, nints, ngroups + + # Check that mask matches 2D data shape + if background_mask.shape[-2:] != image_shape: + log.warning("Mask does not match data shape. Step will be skipped.") + return mismatch, ndim, nints, ngroups + + return False, ndim, nints, ngroups + + +def _clean_one_image(image, mask, background_method, background_box_size, + n_sigma, fit_method, detector, fc, axis_to_correct, + fit_by_channel): + """ + Clean an image by fitting and removing background noise. + + Parameters + ---------- + image : array-like of float + The input image. + mask : array-like of bool + The scene mask. All non-background signal should be + marked as True. + background_method : str + The method for fitting the background. + background_box_size : tuple of int + Background box size, used with `background_method` + is 'model'. + n_sigma : float + N-sigma rejection level for outliers. + fit_method : str + The method for fitting the noise after background + signal is removed. + detector : str + The detector name. + fc : tuple of int + Frequency cutoff values, used for `fit_method` = 'fft'. + axis_to_correct : int + Axis along which noise appears, used for + `fit_method` = 'median'. + fit_by_channel : bool + Option to fit noise in detector channels separately, + used for `fit_method` = 'median' + + Returns + ------- + cleaned_image : array-like of float + The cleaned image. + background : array-like of float + The background fit and removed prior to cleaning. + Used for diagnostic purposes. + success : bool + True if cleaning proceeded as expected; False if + cleaning failed and the step should be skipped. + """ + success = True + + # Make sure data is float32 + image = image.astype(np.float32) + + # Find and replace/mask NaNs + nan_pix = np.isnan(image) + image[nan_pix] = 0.0 + mask[nan_pix] = False + + # If there's no good data remaining in the mask, + # skip this image + if not np.any(mask): + return None, None, success + + # Fit and remove a background level + if str(background_method).lower() == 'none': + background = 0.0 + bkg_sub = image + else: + background = background_level( + image, mask, background_method=background_method, + background_box_size=background_box_size) + log.debug(f'Background level: {np.nanmedian(background):.5g}') + bkg_sub = image - background + + # Flag more signal in the background subtracted image, + # with sigma set by the lower half of the distribution only + clip_to_background( + bkg_sub, mask, sigma_lower=n_sigma, + sigma_upper=n_sigma, lower_half_only=True) + + # Clean the noise + if fit_method == 'fft': + if bkg_sub.shape == (2048, 2048): + # Clean a full-frame image + cleaned_image = fft_clean_full_frame(bkg_sub, mask, detector) + else: + # Clean a subarray image + cleaned_image = fft_clean_subarray(bkg_sub, mask, detector, fc=fc) + else: + cleaned_image = median_clean(bkg_sub, mask, axis_to_correct, + fit_by_channel=fit_by_channel) + if cleaned_image is None: + return None, None, False + + # Restore the background level + cleaned_image += background + + # Restore NaNs in cleaned image + cleaned_image[nan_pix] = np.nan + + return cleaned_image, background, success + + +def _mask_unusable(mask, dq): + """ + Mask unusable data, according to DQ. + + Currently, JUMP and DO_NOT_USE flags are used. + The mask is updated in place + + Parameters + ---------- + mask : array-like of bool + Input mask, updated in place. + dq : array-like of int + DQ flag array matching the mask shape. + """ + dnu = (dq & dqflags.group['DO_NOT_USE']) > 0 + mask[dnu] = False + jump = (dq & dqflags.group['JUMP_DET']) > 0 + mask[jump] = False + + +def do_correction(input_model, input_dir=None, fit_method='median', + fit_by_channel=False, background_method='median', + background_box_size=None, + mask_science_regions=False, n_sigma=2.0, fit_histogram=False, + single_mask=True, user_mask=None, save_mask=False, + save_background=False, save_noise=False): + """ + Apply the 1/f noise correction. + + Parameters + ---------- + input_model : `~jwst.datamodel.JwstDataModel` + Science data to be corrected. + + input_dir : str + Path to the input directory. Used by sub-steps (e.g. assign_wcs + for NIRSpec MOS data) to find auxiliary data. + + fit_method : {'fft', 'median'}, optional + The algorithm to use to fit background noise. + + fit_by_channel : bool, optional + If set, flicker noise is fit independently for each detector channel. + Ignored for MIRI, for subarray data, and for `fit_method` = 'fft'. + + background_method : {'median', 'model', None}, optional + If 'median', the preliminary background to remove and restore + is a simple median of the background data. If 'model', the + background data is fit with a low-resolution model via + `~photutils.background.Background2D`. If None, the background + value is 0.0. + + background_box_size : tuple of int, optional + Box size for the data grid used by `Background2D` when + `background_method` = 'model'. For best results, use a box size + that evenly divides the input image shape. + + mask_science_regions : bool, optional + For NIRSpec, mask regions of the image defined by WCS bounding + boxes for slits/slices, as well as any regions known to be + affected by failed-open MSA shutters. For MIRI imaging, mask + regions of the detector not used for science. + + n_sigma : float, optional + N-sigma rejection level for finding outliers. + + fit_histogram : bool, optional + If set, the 'sigma' used with `n_sigma` for clipping outliers + is derived from a Gaussian fit to a histogram of values. + Otherwise, a simple iterative sigma clipping is performed. + + single_mask : bool, optional + If set, a single mask will be created, regardless of + the number of input integrations. Otherwise, the mask will + be a 3D cube, with one plane for each integration. + + user_mask : str or None, optional + Path to user-supplied mask image. + + save_mask : bool, optional + Switch to indicate whether the mask should be saved. + + save_background : bool, optional + Switch to indicate whether the fit background should be saved. + + save_noise : bool, optional + Switch to indicate whether the fit noise should be saved. + + Returns + ------- + output_model : `~jwst.datamodel.JwstDataModel` + Corrected data. + + mask_model : `~jwst.datamodel.JwstDataModel` + Pixel mask to be saved or None. + + background_model : `~jwst.datamodel.JwstDataModel` + Background model to be saved or None. + + noise_model : `~jwst.datamodel.JwstDataModel` + Background model to be saved or None. + + status : {'COMPLETE', 'SKIPPED'} + Completion status. If errors were encountered, status = 'SKIPPED' + and the output data matches the input data. Otherwise, + status = 'COMPLETE'. + """ + # Track the completion status, for various failure conditions + status = 'SKIPPED' + + detector = input_model.meta.instrument.detector.upper() + subarray = input_model.meta.subarray.name.upper() + exp_type = input_model.meta.exposure.type + slowaxis = input_model.meta.subarray.slowaxis + log.info(f'Input exposure type is {exp_type}, detector={detector}') + + # Check for a valid input that we can work on + if not _check_input(exp_type, fit_method): + return input_model, None, None, None, status + + output_model = input_model.copy() + + # Get parameters needed for subsequent corrections, as appropriate + # to the input data + axis_to_correct, background_method, fit_by_channel, fc = _standardize_parameters( + exp_type, subarray, slowaxis, background_method, fit_by_channel) + + # Make a rate file if needed + if user_mask is None: + image_model = _make_processed_rate_image( + input_model, single_mask, input_dir, exp_type, mask_science_regions) + else: + image_model = input_model + + # Make a mask model from the user input or the rate data + background_mask, mask_model = _make_scene_mask( + user_mask, image_model, mask_science_regions, + n_sigma, fit_histogram, single_mask, save_mask) + + log.info(f"Cleaning image {input_model.meta.filename}") + + # Check data shapes for 2D, 3D, or 4D inputs + mismatch, ndim, nints, ngroups = _check_data_shapes( + input_model, background_mask) + if mismatch: + return output_model, None, None, None, status + + # Close the draft rate model if created - it is no longer needed. + if image_model is not input_model: + image_model.close() + del image_model + + # Make a background cube for saving, if desired + if save_background: + background_to_save = np.zeros_like(input_model.data) + else: + background_to_save = None + + # Loop over integrations and groups (even if there's only 1) + for i in range(nints): + log.debug(f"Working on integration {i + 1}") + for j in range(ngroups): + log.debug(f"Working on group {j + 1}") + + # Copy the scene mask, for further flagging + if background_mask.ndim == 3: + mask = background_mask[i].copy() + else: + mask = background_mask.copy() + + # Get the relevant image data + if ndim == 2: + image = input_model.data + elif ndim == 3: + image = input_model.data[i] + else: + # Ramp data input: + # subtract the current group from the next one + image = input_model.data[i, j+1] - input_model.data[i, j] + dq = input_model.groupdq[i, j+1] + + # Mask any DNU and JUMP pixels + _mask_unusable(mask, dq) + + # Clean the image + cleaned_image, background, success = _clean_one_image( + image, mask, background_method, background_box_size, n_sigma, + fit_method, detector, fc, axis_to_correct, fit_by_channel) + + if not success: + # Cleaning failed for internal reasons - probably the + # mask is not a good match to the data. + log.error(f'Cleaning failed for integration {i + 1}, group {j + 1}') + + # Restore input data to make sure any partial changes + # are thrown away + output_model.data = input_model.data.copy() + return output_model, None, None, None, status + + if cleaned_image is None: + # Cleaning did not proceed because the image is bad: + # leave it as is but continue correcting the rest + log.warning(f'No usable data in integration {i+1}, group {j+1}. ' + f'Skipping correction for this image.') + continue + + # Store the cleaned image in the output model + if ndim == 2: + output_model.data = cleaned_image + if save_background: + background_to_save[:] = background + elif ndim == 3: + output_model.data[i] = cleaned_image + if save_background: + background_to_save[i] = background + else: + # Add the cleaned data diff to the previously cleaned group, + # rather than the noisy input group + output_model.data[i, j+1] = output_model.data[i, j] + cleaned_image + if save_background: + background_to_save[i, j+1] = background + + # Store the background image in a model, if requested + if save_background: + background_model = _make_intermediate_model(output_model, background_to_save) + else: + background_model = None + + # Make a fit noise model for diagnostic purposes by + # diffing the input and output models + if save_noise: + noise_data = output_model.data - input_model.data + noise_model = _make_intermediate_model(output_model, noise_data) + else: + noise_model = None + + # Set completion status + status = 'COMPLETE' + + return output_model, mask_model, background_model, noise_model, status diff --git a/jwst/clean_flicker_noise/clean_flicker_noise_step.py b/jwst/clean_flicker_noise/clean_flicker_noise_step.py new file mode 100644 index 0000000000..a61c179bc8 --- /dev/null +++ b/jwst/clean_flicker_noise/clean_flicker_noise_step.py @@ -0,0 +1,146 @@ +from stdatamodels.jwst import datamodels + +from ..stpipe import Step +from . import clean_flicker_noise + +__all__ = ["CleanFlickerNoiseStep"] + + +class CleanFlickerNoiseStep(Step): + """ + CleanFlickerNoiseStep: This step performs flicker noise correction ("cleaning"). + + Input data is expected to be a ramp file (RampModel), in between + jump and ramp fitting steps, or a rate file (ImageModel or CubeModel). + + Correction algorithms implemented are: + - `fft`: Background noise is fit in frequency space. + Implementation is based on the NSClean algorithm, developed + by Bernard Rauscher. + - `median`: Background noise is characterized by a median + along the detector slow axis. Implementation is based on the + `image1overf` algorithm, developed by Chris Willott. + """ + + class_alias = "clean_flicker_noise" + + spec = """ + fit_method = option('fft', 'median', default='median') # Noise fitting algorithm + fit_by_channel = boolean(default=False) # Fit noise separately by amplifier (NIR only) + background_method = option('median', 'model', None, default='median') # Background fitting algorithm + background_box_size = int_list(min=2, max=2, default=None) # Background box size for modeled background + mask_science_regions = boolean(default=False) # Mask known science regions + n_sigma = float(default=2.0) # Clipping level for non-background signal + fit_histogram = boolean(default=False) # Fit a value histogram to derive sigma + single_mask = boolean(default=True) # Make a single mask for all integrations + user_mask = string(default=None) # Path to user-supplied mask + save_mask = boolean(default=False) # Save the created mask + save_background = boolean(default=False) # Save the fit background + save_noise = boolean(default=False) # Save the fit noise + skip = boolean(default=True) # By default, skip the step + """ + + def process(self, input): + """ + Fit and subtract 1/f background noise from a ramp data set. + + Parameters + ---------- + input : DataModel + Input datamodel to be corrected + + fit_method : str, optional + The noise fitting algorithm to use. Options are 'fft' and 'median'. + + fit_by_channel : bool, optional + If set, flicker noise is fit independently for each detector channel. + Ignored for MIRI, for subarray data, and for `fit_method` = 'fft' + + background_method : {'median', 'model', None} + If 'median', the preliminary background to remove and restore + is a simple median of the background data. If 'model', the + background data is fit with a low-resolution model via + `~photutils.background.Background2D`. If None, the background + value is 0.0. + + background_box_size : tuple of int, optional + Box size for the data grid used by `Background2D` when + `background_method` = 'model'. For best results, use a box size + that evenly divides the input image shape. + + mask_science_regions : bool, optional + For NIRSpec, mask regions of the image defined by WCS bounding + boxes for slits/slices, as well as any regions known to be + affected by failed-open MSA shutters. For MIRI imaging, mask + regions of the detector not used for science. + + n_sigma : float, optional + Sigma clipping threshold to be used in detecting outliers in the image. + + fit_histogram : bool, optional + If set, the 'sigma' used with `n_sigma` for clipping outliers + is derived from a Gaussian fit to a histogram of values. + Otherwise, a simple iterative sigma clipping is performed. + + single_mask : bool, optional + If set, a single mask will be created, regardless of + the number of input integrations. Otherwise, the mask will + be a 3D cube, with one plane for each integration. + + user_mask : None, str, or `~jwst.datamodels.ImageModel` + Optional user-supplied mask image; path to file or opened datamodel. + + save_mask : bool, optional + Save the computed mask image. + + save_background : bool, optional + Save the computed background image. + + save_noise : bool, optional + Save the computed noise image. + + Returns + ------- + output_model : DataModel + The flicker noise corrected datamodel + """ + + # Open the input data model + with datamodels.open(input) as input_model: + + result = clean_flicker_noise.do_correction( + input_model, self.input_dir, self.fit_method, self.fit_by_channel, + self.background_method, self.background_box_size, + self.mask_science_regions, self.n_sigma, self.fit_histogram, + self.single_mask, self.user_mask, + self.save_mask, self.save_background, self.save_noise) + output_model, mask_model, background_model, noise_model, status = result + + # Save the mask, if requested + if self.save_mask and mask_model is not None: + mask_path = self.make_output_path( + basepath=input_model.meta.filename, suffix='mask') + self.log.info(f"Saving mask file {mask_path}") + mask_model.save(mask_path) + mask_model.close() + + # Save the background, if requested + if self.save_background and background_model is not None: + bg_path = self.make_output_path( + basepath=input_model.meta.filename, suffix='flicker_bkg') + self.log.info(f"Saving background file {bg_path}") + background_model.save(bg_path) + background_model.close() + + # Save the noise, if requested + if self.save_noise and noise_model is not None: + noise_path = self.make_output_path( + basepath=input_model.meta.filename, suffix='flicker_noise') + self.log.info(f"Saving noise file {noise_path}") + noise_model.save(noise_path) + noise_model.close() + + # Set the step completion status + output_model.meta.cal_step.clean_flicker_noise = status + + return output_model diff --git a/jwst/nsclean/lib.py b/jwst/clean_flicker_noise/lib.py similarity index 99% rename from jwst/nsclean/lib.py rename to jwst/clean_flicker_noise/lib.py index d4bc05c0fd..0af311dd26 100644 --- a/jwst/nsclean/lib.py +++ b/jwst/clean_flicker_noise/lib.py @@ -517,7 +517,6 @@ def fit(self, return_fit=False, weight_fit=False): if return_fit: return(rfft) - def clean(self, weight_fit=True, return_model=False): """ Clean the data @@ -533,7 +532,7 @@ def clean(self, weight_fit=True, return_model=False): Returns ------- - data : 2D float array + data : array-like of float The cleaned data array. """ self.fit(weight_fit=weight_fit) # Fit the background model diff --git a/jwst/clean_flicker_noise/tests/__init__.py b/jwst/clean_flicker_noise/tests/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/jwst/clean_flicker_noise/tests/test_clean_flicker_noise.py b/jwst/clean_flicker_noise/tests/test_clean_flicker_noise.py new file mode 100644 index 0000000000..7bc237bc94 --- /dev/null +++ b/jwst/clean_flicker_noise/tests/test_clean_flicker_noise.py @@ -0,0 +1,844 @@ +import logging + +import gwcs +import numpy as np +import pytest +from stdatamodels.jwst import datamodels + +from jwst.assign_wcs.tests.test_nirspec import ( + create_nirspec_ifu_file, create_nirspec_fs_file) +from jwst.msaflagopen.tests.test_msa_open import make_nirspec_mos_model, get_file_path +from jwst.clean_flicker_noise import clean_flicker_noise as cfn +from jwst.tests.helpers import LogWatcher + + +def add_metadata(model, shape): + model.meta.instrument.name = 'MIRI' + model.meta.instrument.detector = 'MIRIMAGE' + model.meta.instrument.filter = 'F480M' + model.meta.observation.date = '2015-10-13' + model.meta.observation.time = "00:00:00" + model.meta.exposure.type = 'MIR_IMAGE' + model.meta.exposure.group_time = 1.0 + model.meta.subarray.name = 'FULL' + model.meta.subarray.xstart = 1 + model.meta.subarray.ystart = 1 + model.meta.subarray.xsize = shape[3] + model.meta.subarray.ysize = shape[2] + model.meta.exposure.frame_time = 1.0 + model.meta.exposure.ngroups = shape[1] + model.meta.exposure.group_time = 1.0 + model.meta.exposure.nints = shape[0] + model.meta.exposure.nframes = 1 + model.meta.exposure.groupgap = 0 + model.meta.exposure.readpatt = 'FASTR1' + model.meta.subarray.slowaxis = 2 + + +def make_small_ramp_model(shape=(3, 5, 10, 10)): + rampmodel = datamodels.RampModel(shape) + add_metadata(rampmodel, shape) + + # Make data with a constant rate + for group in range(shape[1]): + rampmodel.data[:, group, :, :] = group + + return rampmodel + + +def make_small_rate_model(shape=(3, 5, 10, 10)): + ratemodel = datamodels.ImageModel(shape[2:]) + add_metadata(ratemodel, shape) + ratemodel.data[:] = 1.0 + return ratemodel + + +def make_small_rateints_model(shape=(3, 5, 10, 10)): + ratemodel = datamodels.CubeModel((shape[0], shape[2], shape[3])) + add_metadata(ratemodel, shape) + ratemodel.data[:] = 1.0 + return ratemodel + + +def make_nirspec_ifu_model(shape=(2048, 2048)): + hdul = create_nirspec_ifu_file(grating='PRISM', filter='CLEAR', + gwa_xtil=0.35986012, gwa_ytil=0.13448857, + gwa_tilt=37.1) + hdul['SCI'].data = np.ones(shape, dtype=float) + rate_model = datamodels.IFUImageModel(hdul) + hdul.close() + return rate_model + + +def make_nirspec_mos_fs_model(): + mos_model = make_nirspec_mos_model() + mos_model.meta.instrument.msa_metadata_file = get_file_path( + 'msa_fs_configuration.fits') + return mos_model + + +def make_nirspec_fs_model(): + hdul = create_nirspec_fs_file(grating="G140M", filter="F100LP") + hdul['SCI'].data = np.ones((2048, 2048), dtype=float) + rate_model = datamodels.ImageModel(hdul) + hdul.close() + + # add the slow axis + rate_model.meta.subarray.slowaxis = 1 + return rate_model + + +class MockUpdate: + def __init__(self): + self.seen = False + + def __call__(self, input_model, mask): + self.seen = True + return mask + + +@pytest.fixture +def log_watcher(monkeypatch): + # Set a log watcher to check for a log message at any level + # in the clean_flicker_noise module + watcher = LogWatcher('') + logger = logging.getLogger('jwst.clean_flicker_noise.clean_flicker_noise') + for level in ['debug', 'info', 'warning', 'error']: + monkeypatch.setattr(logger, level, watcher) + return watcher + + +def test_make_rate(log_watcher): + shape = (3, 5, 10, 10) + ramp_model = make_small_ramp_model(shape) + + log_watcher.message = 'Creating draft rate file' + result = cfn.make_rate(ramp_model, return_cube=True) + assert isinstance(result, datamodels.CubeModel) + assert result.data.shape == (shape[0], shape[2], shape[3]) + assert np.all(result.data == 1.0) + + result = cfn.make_rate(ramp_model, return_cube=False) + assert isinstance(result, datamodels.ImageModel) + assert result.data.shape == shape[2:] + assert np.all(result.data == 1.0) + + ramp_model.close() + result.close() + + # Check for expected log message + log_watcher.assert_seen() + + +def test_postprocess_rate_nirspec(log_watcher): + rate_model = make_nirspec_mos_model() + + log_watcher.message = 'Assigning a WCS' + result = cfn.post_process_rate(rate_model, assign_wcs=True) + assert isinstance(result.meta.wcs, gwcs.WCS) + assert np.sum(result.dq & datamodels.dqflags.pixel['MSA_FAILED_OPEN']) == 0 + log_watcher.assert_seen() + + log_watcher.message = 'Flagging failed-open' + result = cfn.post_process_rate(result, msaflagopen=True) + assert np.sum(result.dq & datamodels.dqflags.pixel['MSA_FAILED_OPEN']) > 0 + log_watcher.assert_seen() + + rate_model.close() + result.close() + + +def test_postprocess_rate_miri(log_watcher): + rate_model = make_small_rate_model() + assert np.sum(rate_model.dq & datamodels.dqflags.pixel['DO_NOT_USE']) == 0 + + log_watcher.message = 'Retrieving flat DQ' + result = cfn.post_process_rate(rate_model, flat_dq=True) + log_watcher.assert_seen() + assert np.sum(result.dq & datamodels.dqflags.pixel['DO_NOT_USE']) > 0 + assert np.all(result.data == rate_model.data) + + rate_model.close() + result.close() + + +@pytest.mark.slow +def test_mask_ifu_slices(): + rate_model = make_nirspec_ifu_model() + + rate_model = cfn.post_process_rate(rate_model, assign_wcs=True) + mask = np.full_like(rate_model.data, True) + + # Mark IFU science regions as False: about 10% of the array + cfn.mask_ifu_slices(rate_model, mask) + assert np.allclose(np.sum(mask), 0.9 * mask.size, atol=0.01 * mask.size) + + rate_model.close() + + +@pytest.mark.parametrize('exptype,blocked', + [('mos', .012), ('mos_fs', .025), ('fs', .05)]) +def test_mask_slits(exptype, blocked): + if exptype == 'mos': + rate_model = make_nirspec_mos_model() + elif exptype == 'mos_fs': + rate_model = make_nirspec_mos_fs_model() + else: + rate_model = make_nirspec_fs_model() + + # Assign a WCS + rate_model = cfn.post_process_rate(rate_model, assign_wcs=True) + + # Mark slits as False + mask = np.full_like(rate_model.data, True) + cfn.mask_slits(rate_model, mask) + + # Check that the fraction of the array blocked is as expected + assert np.allclose(np.sum(mask) / mask.size, 1 - blocked, atol=0.001) + + rate_model.close() + + +@pytest.mark.parametrize('fit_histogram', [True, False]) +@pytest.mark.parametrize('lower_half_only', [True, False]) +def test_clip_to_background(log_watcher, fit_histogram, lower_half_only): + """Integrity checks for clipping a simple data array.""" + # Make an array with normal noise + shape = (100, 100) + rng = np.random.default_rng(seed=123) + image = rng.normal(0, 0.01, size=shape) + mask = np.full(shape, True) + + # Add a strong high outlier, strong low outlier + image[0, 0] += 1 + image[1, 1] += -1 + + # Center found is close to zero, printed when verbose=True + log_watcher.message = "center: 0.000" + cfn.clip_to_background( + image, mask, fit_histogram=fit_histogram, + lower_half_only=lower_half_only, verbose=True) + log_watcher.assert_seen() + + # All outliers clipped with defaults, as well as a small + # percent of the remaining data + assert not mask[0, 0] + assert not mask[1, 1] + assert np.allclose(np.sum(mask) / mask.size, 0.98, atol=0.019) + + # Upper outlier stays with large sigma_upper + mask = np.full(shape, True) + cfn.clip_to_background( + image, mask, fit_histogram=fit_histogram, + lower_half_only=lower_half_only, sigma_upper=1000) + assert mask[0, 0] + assert not mask[1, 1] + assert np.allclose(np.sum(mask) / mask.size, 0.98, atol=0.019) + + # Lower outlier stays with large sigma_lower + mask = np.full(shape, True) + cfn.clip_to_background( + image, mask, fit_histogram=fit_histogram, + lower_half_only=lower_half_only, sigma_lower=1000) + assert not mask[0, 0] + assert mask[1, 1] + assert np.allclose(np.sum(mask) / mask.size, 0.98, atol=0.019) + + +def test_clip_to_background_fit_fails(log_watcher): + shape = (10, 10) + + # histogram failure: all data NaN + log_watcher.message = "Histogram failed" + image = np.full(shape, np.nan) + mask = np.full(shape, True) + with pytest.warns(RuntimeWarning): + cfn.clip_to_background(image, mask, fit_histogram=True, verbose=True) + assert np.all(mask) + log_watcher.assert_seen() + + # if mask is all False, warning is avoided, mask is unchanged + mask[:] = False + cfn.clip_to_background(image, mask, fit_histogram=True, verbose=True) + assert not np.all(mask) + + # fit failure: data is not normal + log_watcher.message = "Gaussian fit failed" + image = np.full(shape, 0.0) + mask = np.full(shape, True) + image[5:, 5:] = 0.1 + cfn.clip_to_background(image, mask, fit_histogram=True, verbose=True) + assert np.all(mask) + log_watcher.assert_seen() + + +@pytest.mark.parametrize('exptype', ['mos', 'mos_fs', 'ifu']) +def test_create_mask_nirspec(monkeypatch, exptype): + # monkeypatch local functions for speed and check that they are called + # (actual behavior tested in separate unit tests) + mock = MockUpdate() + if exptype == 'mos': + # NIRSpec MOS data + monkeypatch.setattr(cfn, 'mask_slits', mock) + rate_model = make_nirspec_mos_model() + elif exptype == 'mos_fs': + monkeypatch.setattr(cfn, 'mask_slits', mock) + rate_model = make_nirspec_mos_fs_model() + + # also assign a WCS so the FS can be detected + rate_model = cfn.post_process_rate(rate_model, assign_wcs=True) + else: + # NIRSpec IFU data + monkeypatch.setattr(cfn, 'mask_ifu_slices', mock) + rate_model = make_nirspec_ifu_model((2048, 2038)) + + rate_model.data[:] = 1.0 + mask = cfn.create_mask(rate_model) + assert mask.shape == rate_model.data.shape + + # Nothing to mask in uniform data + assert np.all(mask) + + # Slit function not called if mask_science_regions not set + assert not mock.seen + + # Fixed slit region not blocked + assert np.all(mask[cfn.NRS_FS_REGION[0]:cfn.NRS_FS_REGION[1]]) + + # Call again but block science regions + mask = cfn.create_mask(rate_model, mask_science_regions=True) + assert mock.seen + if exptype == 'mos_fs': + # FS region still not blocked - may need correction + assert np.all(mask[cfn.NRS_FS_REGION[0]:cfn.NRS_FS_REGION[1]]) + else: + assert not np.any(mask[cfn.NRS_FS_REGION[0]:cfn.NRS_FS_REGION[1]]) + + +def test_create_mask_miri(): + # small MIRI imaging data + rate_model = make_small_rate_model() + rate_model.dq[5, 5] = datamodels.dqflags.pixel['DO_NOT_USE'] + + mask = cfn.create_mask(rate_model) + assert mask.shape == rate_model.data.shape + + # Nothing to mask in uniform data + assert np.all(mask) + + # Call again but block non-science regions + mask = cfn.create_mask(rate_model, mask_science_regions=True) + assert not mask[5, 5] + assert np.sum(mask) == mask.size - 1 + + rate_model.close() + + +def test_create_mask_from_rateints(): + # small rateints data + rate_model = make_small_rateints_model() + + # Add an outlier in each integration + for i in range(rate_model.data.shape[0]): + rate_model.data[i, i, i] += 100 + + mask = cfn.create_mask(rate_model) + assert mask.shape == rate_model.data.shape + assert mask.ndim == 3 + + # Outlier is masked in each integration + assert not mask[0, 0, 0] + assert not mask[1, 1, 1] + assert not mask[2, 2, 2] + assert np.sum(mask) == mask.size - 3 + + # Call again but make a single mask: + # bad data mask is or-ed across integrations + mask = cfn.create_mask(rate_model, single_mask=True) + assert mask.shape == rate_model.data.shape[1:] + assert mask.ndim == 2 + + assert not mask[0, 0] + assert not mask[1, 1] + assert not mask[2, 2] + assert np.sum(mask) == mask.size - 3 + + rate_model.close() + + +def test_background_level(log_watcher): + shape = (100, 100) + image = np.full(shape, 1.0) + mask = np.full(shape, True) + + # add an outlier to be clipped + image[50, 50] = 1000 + + # no background + background = cfn.background_level(image, mask, background_method=None) + assert background == 0.0 + + # median method + background = cfn.background_level(image, mask, background_method='median') + assert background == 1.0 + + # model method + background = cfn.background_level( + image, mask, background_method='model', background_box_size=(10, 10)) + assert background.shape == shape + assert np.all(background == 1.0) + + # model method with mismatched box size: + # warns, but completes successfully + log_watcher.message = "does not divide evenly" + background = cfn.background_level( + image, mask, background_method='model', background_box_size=None) + assert background.shape == shape + assert np.all(background == 1.0) + log_watcher.assert_seen() + + # make image mostly bad, one good region + image[:] = np.nan + image[20:25, 20:25] = 1.0 + + # background fit fails: falls back on simple median + log_watcher.message = "Background fit failed, using median" + background = cfn.background_level( + image, mask, background_method='model', background_box_size=(10, 10)) + assert background == 1.0 + log_watcher.assert_seen() + + +@pytest.mark.parametrize('array_type,fraction_good', [('full', 0.65), ('subarray', 0.37)]) +@pytest.mark.parametrize('detector', ['NRS1', 'NRS2']) +def test_fft_clean(array_type, fraction_good, detector): + if array_type == 'full': + clean_function = cfn.fft_clean_full_frame + else: + clean_function = cfn.fft_clean_subarray + + shape = (80, 80) + mask = np.full(shape, True) + + # zero image should still come out zero + image = np.full(shape, 0.0) + cleaned_image = clean_function(image.copy(), mask, detector) + assert np.allclose(cleaned_image, 0.0) + + # image with regular vertical pattern, centered on 0 + high = np.full((80, 16), 0.1) + low = np.full((80, 16), -0.1) + image = np.hstack([high, low, high, low, high]) + cleaned_image = clean_function(image.copy(), mask, detector) + + # results should be mostly close to zero, + # some artifacts at pattern and edge boundaries, + # edge boundaries are currently worse for subarray correction + good_correction = np.abs(cleaned_image) < 0.01 + assert np.allclose(np.sum(good_correction) / cleaned_image.size, + fraction_good, atol=0.1) + + +def test_fft_clean_error(monkeypatch, log_watcher): + def raise_error(*args, **kwargs): + raise np.linalg.LinAlgError('Linear algebra error') + + clean_function = cfn.fft_clean_full_frame + monkeypatch.setattr(cfn.NSClean, 'clean', raise_error) + + shape = (10, 10) + mask = np.full(shape, True) + image = np.full(shape, 1.0) + + log_watcher.message = "Error cleaning image" + cleaned_image = clean_function(image.copy(), mask, 'NRS1') + assert cleaned_image is None + log_watcher.assert_seen() + + +def test_fft_subarray_clean_error(monkeypatch, log_watcher): + shape = (10, 10) + image = np.full(shape, 1.0) + + # Mask is all bad: error message, returns None + mask = np.full(shape, False) + log_watcher.message = "No good pixels" + cleaned_image = cfn.fft_clean_subarray(image.copy(), mask, 'NRS1') + assert cleaned_image is None + log_watcher.assert_seen() + + # Mask is nearly all bad except a pixel at the edge: + # triggers a linear algebra error + mask[0, 0] = True + log_watcher.message = "Error cleaning image" + cleaned_image = cfn.fft_clean_subarray(image.copy(), mask, 'NRS1') + assert cleaned_image is None + log_watcher.assert_seen() + + # Mask is mostly bad: warns but continues + mask[5, 5] = True + log_watcher.message = "Insufficient reference pixels" + cleaned_image = cfn.fft_clean_subarray(image.copy(), mask, 'NRS1', minfrac=0.5) + assert np.allclose(cleaned_image, image) + log_watcher.assert_seen() + + +def test_median_clean(): + shape = (20, 20) + mask = np.full(shape, True) + + # zero image should still come out zero + image = np.full(shape, 0.0) + cleaned_image = cfn.median_clean(image, mask, 1) + assert np.allclose(cleaned_image, 0.0) + + # image with regular vertical pattern, centered on 0 + high = np.full((20, 4), 0.1) + low = np.full((20, 4), -0.1) + image = np.hstack([high, low, high, low, high]) + + # clean along vertical axis - + # image should be all zero + cleaned_image = cfn.median_clean(image, mask, 1) + assert np.allclose(cleaned_image, 0.0) + + # clean along horizontal axis - + # cleaning removes median value, from high stripes + cleaned_image = cfn.median_clean(image, mask, 2) + assert np.allclose(cleaned_image, image - 0.1) + + +def test_median_clean_by_channel(): + shape = (2048, 2048) + mask = np.full(shape, True) + + # zero image should still come out zero + image = np.full(shape, 0.0) + cleaned_image = cfn.median_clean(image, mask, 1, fit_by_channel=True) + assert np.allclose(cleaned_image, 0.0) + + # image with regular vertical pattern, centered on 0 + high = np.full((2048, 16), 0.1) + low = np.full((2048, 16), -0.1) + image = np.hstack([high, low] * 64) + + # add offset by vertical channel + offset = np.zeros_like(image) + offset[:512] += 0.2 + offset[512:1024] += 0.3 + offset[1024:1536] += 0.4 + offset[1536:] += 0.5 + image += offset + + # clean along vertical axis by channel - + # image should be all zero + cleaned_image = cfn.median_clean(image, mask, 1, fit_by_channel=True) + assert np.allclose(cleaned_image, 0.0) + + # clean along vertical axis for the whole image - + # vertical pattern is removed, offset by channel stays, + # but median value is subtracted + cleaned_image = cfn.median_clean(image, mask, 1, fit_by_channel=False) + assert np.allclose(cleaned_image, offset - 0.35) + + # clean along horizontal axis - + # cleaning removes offset value, leaves vertical pattern + cleaned_image = cfn.median_clean(image, mask, 2, fit_by_channel=True) + assert np.allclose(cleaned_image, image - offset) + + +@pytest.mark.parametrize('mask_science', [True, False]) +def test_do_correction_miri_imaging_ramp(mask_science): + ramp_model = make_small_ramp_model() + cleaned, mask, bg, noise, status = cfn.do_correction( + ramp_model, mask_science_regions=mask_science) + + # uniform data, correction has no effect + assert cleaned.data.shape == ramp_model.shape + assert np.allclose(cleaned.data, ramp_model.data) + assert status == 'COMPLETE' + + # extra models are not created + assert mask is None + assert bg is None + assert noise is None + + ramp_model.close() + cleaned.close() + + +@pytest.mark.parametrize('mask_science', [True, False]) +def test_do_correction_nirspec_rate(mask_science): + rate_model = make_nirspec_fs_model() + cleaned, mask, bg, noise, status = cfn.do_correction( + rate_model, mask_science_regions=mask_science) + + # uniform data, correction has no effect + assert cleaned.data.shape == rate_model.data.shape + assert np.allclose(cleaned.data, rate_model.data) + assert status == 'COMPLETE' + + # extra models are not created + assert mask is None + assert bg is None + assert noise is None + + rate_model.close() + cleaned.close() + + +@pytest.mark.parametrize('single_mask', [True, False]) +def test_do_correction_rateints(single_mask): + rate_model = make_small_rateints_model() + cleaned, mask, bg, noise, status = cfn.do_correction( + rate_model, single_mask=single_mask, save_mask=True) + + # uniform data, correction has no effect + assert cleaned.data.shape == rate_model.shape + assert np.allclose(cleaned.data, rate_model.data) + assert status == 'COMPLETE' + + if single_mask: + assert mask.shape == cleaned.data.shape[-2:] + else: + assert mask.shape == cleaned.data.shape + + rate_model.close() + cleaned.close() + mask.close() + + +def test_do_correction_unsupported(log_watcher): + ramp_model = make_small_ramp_model() + ramp_model.meta.exposure.type = 'MIR_MRS' + + log_watcher.message = "not supported" + cleaned, _, _, _, status = cfn.do_correction(ramp_model) + assert cleaned is ramp_model + assert status == 'SKIPPED' + log_watcher.assert_seen() + + ramp_model.close() + + +def test_do_correction_no_fit_by_channel(log_watcher): + ramp_model = make_small_ramp_model() + + # fit_by_channel is only used for NIR data - + # log a warning, but step still completes + log_watcher.message = "can only be used for full-frame NIR" + cleaned, _, _, _, status = cfn.do_correction(ramp_model, fit_by_channel=True) + assert status == 'COMPLETE' + log_watcher.assert_seen() + + ramp_model.close() + cleaned.close() + + +@pytest.mark.parametrize('exptype', ['MIR_IMAGE', 'NRC_IMAGE', 'NIS_IMAGE']) +def test_do_correction_fft_not_allowed(log_watcher, exptype): + ramp_model = make_small_ramp_model() + ramp_model.meta.exposure.type = exptype + + # not allowed for MIRI, NIRCAM, NIRISS + log_watcher.message = f"cannot be applied to exp_type {exptype}" + cleaned, _, _, _, status = cfn.do_correction(ramp_model, fit_method='fft') + assert cleaned is ramp_model + assert status == 'SKIPPED' + log_watcher.assert_seen() + + ramp_model.close() + + +@pytest.mark.parametrize('none_value', [None, 'none', 'None']) +def test_do_correction_no_background(none_value): + ramp_model = make_small_ramp_model() + + # Don't remove background before fitting noise + cleaned, _, _, _, _ = cfn.do_correction( + ramp_model, background_method=none_value) + + # Output data is all zero: uniform level is removed + assert np.allclose(cleaned.data, 0.0) + + ramp_model.close() + cleaned.close() + + +@pytest.mark.parametrize('ndim', [2, 3]) +def test_do_correction_user_mask(tmp_path, ndim): + ramp_model = make_small_ramp_model() + + if ndim == 3: + mask = np.full((ramp_model.shape[0], + ramp_model.shape[2], + ramp_model.shape[3]), True) + mask_model = datamodels.CubeModel(mask) + else: + mask = np.full(ramp_model.shape[-2:], True) + mask_model = datamodels.ImageModel(mask) + user_mask = str(tmp_path / 'mask.fits') + mask_model.save(user_mask) + mask_model.close() + + cleaned, output_mask, _, _, _ = cfn.do_correction( + ramp_model, user_mask=user_mask, save_mask=True) + + assert np.all(output_mask.data == mask_model.data) + assert np.allclose(cleaned.data, ramp_model.data) + + ramp_model.close() + cleaned.close() + output_mask.close() + + +@pytest.mark.parametrize('input_type', ['rate', 'rateints', 'ramp']) +def test_do_correction_user_mask_mismatch(tmp_path, input_type, log_watcher): + shape = (3, 5, 20, 20) + if input_type == 'rate': + model = make_small_rate_model(shape) + elif input_type == 'rateints': + model = make_small_rateints_model(shape) + else: + model = make_small_ramp_model(shape) + + mask = np.full((10, 10), True) + mask_model = datamodels.ImageModel(mask) + user_mask = str(tmp_path / 'mask.fits') + mask_model.save(user_mask) + mask_model.close() + + log_watcher.message = 'Mask does not match' + cleaned, output_mask, _, _, status = cfn.do_correction( + model, user_mask=user_mask, save_mask=True) + + log_watcher.assert_seen() + assert status == 'SKIPPED' + assert output_mask is None + + model.close() + cleaned.close() + + +def test_do_correction_user_mask_mismatch_integ(tmp_path, log_watcher): + shape = (3, 5, 20, 20) + model = make_small_rateints_model(shape) + + mask = np.full((2, 20, 20), True) + mask_model = datamodels.CubeModel(mask) + user_mask = str(tmp_path / 'mask.fits') + mask_model.save(user_mask) + mask_model.close() + + log_watcher.message = 'Mask does not match' + cleaned, output_mask, _, _, status = cfn.do_correction( + model, user_mask=user_mask, save_mask=True) + + log_watcher.assert_seen() + assert status == 'SKIPPED' + assert output_mask is None + + model.close() + cleaned.close() + + +@pytest.mark.parametrize('subarray', ['SUBS200A1', 'ALLSLITS']) +def test_do_correction_fft_subarray(subarray): + ramp_model = make_small_ramp_model() + ramp_model.meta.exposure.type = 'NRS_FIXEDSLIT' + ramp_model.meta.subarray.slowaxis = 1 + ramp_model.meta.subarray.name = subarray + + cleaned, _, _, noise, _ = cfn.do_correction( + ramp_model, fit_method='fft', save_noise=True) + + # uniform data, correction has no effect + assert cleaned.data.shape == ramp_model.data.shape + assert np.allclose(cleaned.data, ramp_model.data) + assert np.allclose(noise.data, 0.0) + + ramp_model.close() + cleaned.close() + noise.close() + + +def test_do_correction_fft_full(): + rate_model = make_nirspec_fs_model() + + cleaned, _, _, noise, _ = cfn.do_correction( + rate_model, fit_method='fft', save_noise=True) + + # uniform data, correction has no effect + assert cleaned.data.shape == rate_model.data.shape + assert np.allclose(cleaned.data, rate_model.data) + assert np.allclose(noise.data, 0.0) + + rate_model.close() + cleaned.close() + noise.close() + + +def test_do_correction_clean_fails(monkeypatch, log_watcher): + ramp_model = make_small_ramp_model() + ramp_model.meta.exposure.type = 'NRS_FIXEDSLIT' + ramp_model.meta.subarray.slowaxis = 1 + + # Add some noise so that input and output are + # expected to be different + rng = np.random.default_rng(seed=123) + ramp_model.data += rng.normal(0, 0.1, size=ramp_model.data.shape) + cleaned, _, _, _, status = cfn.do_correction( + ramp_model, fit_method='fft') + assert not np.allclose(cleaned.data, ramp_model.data) + assert status == 'COMPLETE' + + # Mock a None-value returned from cleaning function + monkeypatch.setattr(cfn, 'fft_clean_subarray', lambda *args, **kwargs: None) + + # Call again + log_watcher.message = "Cleaning failed" + cleaned, _, _, _, status = cfn.do_correction( + ramp_model, fit_method='fft') + + # Error message issued, status is skipped, + # output data is the same as input + log_watcher.assert_seen() + assert status == 'SKIPPED' + assert np.allclose(cleaned.data, ramp_model.data) + + ramp_model.close() + cleaned.close() + + +@pytest.mark.parametrize('save_type', ['noise', 'background']) +@pytest.mark.parametrize('input_type', ['rate', 'rateints', 'ramp']) +def test_do_correction_save_intermediate(save_type, input_type): + shape = (3, 5, 20, 20) + if input_type == 'rate': + model = make_small_rate_model(shape) + elif input_type == 'rateints': + model = make_small_rateints_model(shape) + else: + model = make_small_ramp_model(shape) + + save_bg = (save_type == 'background') + save_noise = (save_type == 'noise') + cleaned, _, background, noise, status = cfn.do_correction( + model, save_background=save_bg, save_noise=save_noise) + + # Output background model always matches input shape + # and datamodel type + if save_bg: + assert background.data.shape == model.data.shape + assert type(background) is type(model) + background.close() + else: + assert background is None + if save_noise: + assert noise.data.shape == model.data.shape + assert type(noise) is type(model) + noise.close() + else: + assert noise is None + + model.close() diff --git a/jwst/clean_flicker_noise/tests/test_clean_flicker_noise_step.py b/jwst/clean_flicker_noise/tests/test_clean_flicker_noise_step.py new file mode 100644 index 0000000000..d8f2e1dcc6 --- /dev/null +++ b/jwst/clean_flicker_noise/tests/test_clean_flicker_noise_step.py @@ -0,0 +1,82 @@ +import os + +import pytest +from stdatamodels.jwst import datamodels + +from jwst.clean_flicker_noise import CleanFlickerNoiseStep +from .test_clean_flicker_noise import make_small_ramp_model + + +@pytest.mark.parametrize('skip', [True, False]) +def test_output_type(skip): + input_model = make_small_ramp_model() + cleaned = CleanFlickerNoiseStep.call(input_model, skip=skip) + + # output is a ramp model either way + assert isinstance(cleaned, datamodels.RampModel) + + # status may be SKIPPED or COMPLETE + if skip: + assert cleaned.meta.cal_step.clean_flicker_noise == 'SKIPPED' + else: + assert cleaned.meta.cal_step.clean_flicker_noise == 'COMPLETE' + + input_model.close() + cleaned.close() + + +def test_save_mask(tmp_path): + input_model = make_small_ramp_model() + input_model.meta.filename = 'test_jump.fits' + + output_dir = str(tmp_path) + cleaned = CleanFlickerNoiseStep.call( + input_model, skip=False, save_mask=True, output_dir=output_dir) + + mask_file = str(tmp_path / 'test_mask.fits') + assert os.path.isfile(mask_file) + + mask_model = datamodels.open(mask_file) + assert isinstance(mask_model, datamodels.ImageModel) + + input_model.close() + cleaned.close() + mask_model.close() + + +def test_save_background(tmp_path): + input_model = make_small_ramp_model() + input_model.meta.filename = 'test_jump.fits' + + output_dir = str(tmp_path) + cleaned = CleanFlickerNoiseStep.call( + input_model, skip=False, save_background=True, output_dir=output_dir) + + bkg_file = str(tmp_path / 'test_flicker_bkg.fits') + assert os.path.isfile(bkg_file) + + bkg_model = datamodels.open(bkg_file) + assert isinstance(bkg_model, datamodels.RampModel) + + input_model.close() + cleaned.close() + bkg_model.close() + + +def test_save_noise(tmp_path): + input_model = make_small_ramp_model() + input_model.meta.filename = 'test_jump.fits' + + output_dir = str(tmp_path) + cleaned = CleanFlickerNoiseStep.call( + input_model, skip=False, save_noise=True, output_dir=output_dir) + + noise_file = str(tmp_path / 'test_flicker_noise.fits') + assert os.path.isfile(noise_file) + + noise_model = datamodels.open(noise_file) + assert isinstance(noise_model, datamodels.RampModel) + + input_model.close() + cleaned.close() + noise_model.close() diff --git a/jwst/lib/suffix.py b/jwst/lib/suffix.py index 3fff82ba9f..a7af86beb7 100644 --- a/jwst/lib/suffix.py +++ b/jwst/lib/suffix.py @@ -159,6 +159,8 @@ 'spec3pipeline', 'nsclean', 'nscleanstep', + 'clean_flicker_noise', + 'cleanflickernoisestep', 'outlierdetectionstep', 'groupscalestep', 'hlspstep', diff --git a/jwst/msaflagopen/msaflag_open.py b/jwst/msaflagopen/msaflag_open.py index b348ff4fb0..39a49e7b7c 100644 --- a/jwst/msaflagopen/msaflag_open.py +++ b/jwst/msaflagopen/msaflag_open.py @@ -84,10 +84,10 @@ def flag(input_datamodel, failed_slitlets, wcs_refnames): science data with DQ flags of affected modified """ - # # Use the machinery in assign_wcs to create a WCS object for the bad shutters pipeline = slitlets_wcs(input_datamodel, wcs_refnames, failed_slitlets) wcs = WCS(pipeline) + # Create output as a copy of the input science data model so we can overwrite # the wcs with the wcs for the failed open shutters # Have to make sure the EXP_TYPE is NRS_MSASPEC so that nrs_wcs_set_input works, @@ -98,31 +98,30 @@ def flag(input_datamodel, failed_slitlets, wcs_refnames): temporary_copy.meta.exposure.type = 'NRS_MSASPEC' dq_array = input_datamodel.dq for slitlet in failed_slitlets: - # # Pick the WCS for this slitlet from the WCS of the exposure thiswcs = nrs_wcs_set_input(temporary_copy, slitlet.name) - # + # Convert the bounding box for this slitlet to a set of indices to use as a slice xmin, xmax, ymin, ymax = boundingbox_to_indices(temporary_copy, thiswcs.bounding_box) - # + # Make a grid of points within the slice y_indices, x_indices = np.mgrid[ymin:ymax, xmin:xmax] - # + # Calculate the arrays of coordinates for each pixel in the slice coordinate_array = thiswcs(x_indices, y_indices) - # + # The coordinate_array is a tuple of arrays, one for each output coordinate # In this case there should be 3 arrays, one each for RA, Dec and Wavelength # For pixels outside the slitlet, the arrays have NaN - # + # Make a subarray from these coordinate arrays by setting pixels that aren't # NaN to FAILDOPENFLAG, the rest to 0 dq_subarray = wcs_to_dq(coordinate_array, FAILEDOPENFLAG) - # + # bitwise-or this subarray with the slice in the original exposure's DQ array dq_array = or_subarray_with_array(dq_array, dq_subarray, xmin, xmax, ymin, ymax) - # + # Set the dq array of the input datamodel to the corrected dq array input_datamodel.dq = dq_array return input_datamodel @@ -148,7 +147,7 @@ def boundingbox_to_indices(data_model, bounding_box): Range of indices of overlap between science data array and bounding box """ - nrows, ncols = data_model.data.shape + nrows, ncols = data_model.data.shape[-2:] ((x1, x2), (y1, y2)) = bounding_box xmin = int(min(x1, x2)) xmin = max(xmin, 0) @@ -243,5 +242,5 @@ def or_subarray_with_array(dq_array, dq_subarray, xmin, xmax, ymin, ymax): Bitwise-or the slice of the dq array with the section corresponding to the failed open shutter """ - dq_array[ymin:ymax, xmin:xmax] = np.bitwise_or(dq_array[ymin:ymax, xmin:xmax], dq_subarray) + dq_array[..., ymin:ymax, xmin:xmax] = np.bitwise_or(dq_array[..., ymin:ymax, xmin:xmax], dq_subarray) return dq_array diff --git a/jwst/msaflagopen/tests/test_msa_open.py b/jwst/msaflagopen/tests/test_msa_open.py index 131c5eba46..c94f7af81b 100644 --- a/jwst/msaflagopen/tests/test_msa_open.py +++ b/jwst/msaflagopen/tests/test_msa_open.py @@ -28,6 +28,55 @@ def get_file_path(filename): return os.path.join(data_path, filename) +def make_nirspec_mos_model(): + im = ImageModel((2048, 2048)) + im.meta.wcsinfo = { + 'dec_ref': -0.00601415671349804, + 'ra_ref': -0.02073605215697509, + 'roll_ref': -0.0, + 'v2_ref': -453.5134, + 'v3_ref': -373.4826, + 'v3yangle': 0.0, + 'vparity': -1} + + im.meta.instrument = { + 'detector': 'NRS1', + 'filter': 'F100LP', + 'grating': 'G140M', + 'name': 'NIRSPEC', + 'gwa_tilt': 37.0610, + 'gwa_xtilt': 0.0001, + 'gwa_ytilt': 0.0001, + 'msa_metadata_id': 12} + + im.meta.observation = { + 'program_number': '1234', + 'date': '2016-09-05', + 'time': '8:59:37'} + + im.meta.exposure = { + 'duration': 11.805952, + 'end_time': 58119.85416, + 'exposure_time': 11.776, + 'frame_time': 0.11776, + 'group_time': 0.11776, + 'groupgap': 0, + 'integration_time': 11.776, + 'nframes': 1, + 'ngroups': 100, + 'nints': 1, + 'nresets_between_ints': 0, + 'nsamples': 1, + 'readpatt': 'NRSRAPID', + 'sample_time': 10.0, + 'start_time': 58119.8333, + 'type': 'NRS_MSASPEC', + 'zero_frame': False} + im.meta.instrument.msa_metadata_file = get_file_path('msa_configuration.fits') + im.meta.dither.position_number = 1 + return im + + def test_or_subarray_with_array(): """test bitwise or with array and subarray.""" @@ -137,52 +186,7 @@ def test_boundingbox_from_indices(): def test_msaflagopen_step(): - im = ImageModel((2048, 2048)) - im.meta.wcsinfo = { - 'dec_ref': -0.00601415671349804, - 'ra_ref': -0.02073605215697509, - 'roll_ref': -0.0, - 'v2_ref': -453.5134, - 'v3_ref': -373.4826, - 'v3yangle': 0.0, - 'vparity': -1} - - im.meta.instrument = { - 'detector': 'NRS1', - 'filter': 'F100LP', - 'grating': 'G140M', - 'name': 'NIRSPEC', - 'gwa_tilt': 37.0610, - 'gwa_xtilt': 0.0001, - 'gwa_ytilt': 0.0001, - 'msa_metadata_id': 12} - - im.meta.observation = { - 'program_number': '1234', - 'date': '2016-09-05', - 'time': '8:59:37'} - - im.meta.exposure = { - 'duration': 11.805952, - 'end_time': 58119.85416, - 'exposure_time': 11.776, - 'frame_time': 0.11776, - 'group_time': 0.11776, - 'groupgap': 0, - 'integration_time': 11.776, - 'nframes': 1, - 'ngroups': 100, - 'nints': 1, - 'nresets_between_ints': 0, - 'nsamples': 1, - 'readpatt': 'NRSRAPID', - 'sample_time': 10.0, - 'start_time': 58119.8333, - 'type': 'NRS_MSASPEC', - 'zero_frame': False} - im.meta.instrument.msa_metadata_file = get_file_path('msa_configuration.fits') - im.meta.dither.position_number = 1 - + im = make_nirspec_mos_model() im = AssignWcsStep.call(im) result = MSAFlagOpenStep.call(im) diff --git a/jwst/nsclean/nsclean.py b/jwst/nsclean/nsclean.py deleted file mode 100644 index 27c6509548..0000000000 --- a/jwst/nsclean/nsclean.py +++ /dev/null @@ -1,508 +0,0 @@ -import logging -import numpy as np -import gwcs -from astropy.stats import sigma_clipped_stats - -from jwst import datamodels -from jwst.assign_wcs import nirspec -from stdatamodels.jwst.datamodels import dqflags - -from jwst.nsclean.lib import NSClean, NSCleanSubarray - -log = logging.getLogger(__name__) -log.setLevel(logging.DEBUG) - - -def mask_ifu_slices(input_model, mask): - - """Find pixels located within IFU slices and flag them in the - mask, so that they do not get used. - - Parameters - ---------- - input_model : data model object - science data - - mask : 2D bool array - input mask that will be updated - - Returns - ------- - mask : 2D bool array - output mask with additional flags for science pixels - """ - - log.info("Finding slice pixels for an IFU image") - - # Initialize global DQ map to all zero (OK to use) - dqmap = np.zeros_like(input_model.dq) - - # Get the wcs objects for all IFU slices - list_of_wcs = nirspec.nrs_ifu_wcs(input_model) - - # Loop over the IFU slices, finding the valid region for each - for (k, ifu_wcs) in enumerate(list_of_wcs): - - # Construct array indexes for pixels in this slice - x, y = gwcs.wcstools.grid_from_bounding_box(ifu_wcs.bounding_box, - step=(1, 1), - center=True) - # Get the world coords for all pixels in this slice; - # all we actually need are wavelengths - coords = ifu_wcs(x, y) - dq = dqmap[y.astype(int), x.astype(int)] - wl = coords[2] - # set non-NaN wavelength locations as do not use (one) - valid = ~np.isnan(wl) - dq = dq[valid] - x = x[valid] - y = y[valid] - dq[:] = 1 - - # Copy DQ for this slice into global DQ map - dqmap[y.astype(int), x.astype(int)] = dq - - # Now set all non-zero locations in the mask to False (do not use) - mask[dqmap==1] = False - - return mask - - -def mask_slits(input_model, mask): - - """Find pixels located within MOS or fixed slit footprints - and flag them in the mask, so that they do not get used. - - Parameters - ---------- - input_model : data model object - science data - - mask : 2D bool array - input mask that will be updated - - Returns - ------- - mask : 2D bool array - output mask with additional flags for slit pixels - """ - - from jwst.extract_2d.nirspec import offset_wcs - - log.info("Finding slit/slitlet pixels") - - # Get the slit-to-msa frame transform from the WCS object - slit2msa = input_model.meta.wcs.get_transform('slit_frame', 'msa_frame') - - # Loop over the slits, marking all the pixels within each bounding - # box as False (do not use) in the mask. - # Note that for 3D masks (TSO mode), all planes will be set to the same value. - for slit in slit2msa.slits: - slit_wcs = nirspec.nrs_wcs_set_input(input_model, slit.name) - xlo, xhi, ylo, yhi = offset_wcs(slit_wcs) - mask[..., ylo:yhi, xlo:xhi] = False - - return mask - - -def create_mask(input_model, mask_spectral_regions, n_sigma): - """Create the pixel mask needed for setting which pixels to use - for measuring 1/f noise. - - Parameters - ---------- - input_model : data model object - science data - - mask_spectral_regions : bool - mask slit/slice regions defined in WCS - - n_sigma : float - sigma threshold for masking outliers - - Returns - ------- - mask : 2D or 3D bool array - image mask - - nan_pix : array - indexes of image locations with NaN values - """ - exptype = input_model.meta.exposure.type.lower() - - # Initialize mask to all True. Subsequent operations will mask - # out pixels that contain signal. - # Note: mask will be 3D for BOTS mode data - mask = np.full(np.shape(input_model.dq), True) - - # Mask any reference pixels - ref_pix = input_model.dq & dqflags.pixel['REFERENCE_PIXEL'] - mask[ref_pix > 0] = False - - # If IFU, mask all pixels contained in the IFU slices - if exptype == 'nrs_ifu' and mask_spectral_regions: - mask = mask_ifu_slices(input_model, mask) - - # If MOS or FS, mask all pixels affected by open slitlets - if exptype in ['nrs_fixedslit', 'nrs_brightobj', 'nrs_msaspec'] and mask_spectral_regions: - mask = mask_slits(input_model, mask) - - # If IFU or MOS, mask pixels affected by failed-open shutters - if exptype in ['nrs_ifu', 'nrs_msaspec']: - open_pix = input_model.dq & dqflags.pixel['MSA_FAILED_OPEN'] - mask[open_pix > 0] = False - - # Temporarily reset NaN pixels and mask them. - # Save the list of NaN pixel coords, so that they can be reset at the end. - nan_pix = np.isnan(input_model.data) - input_model.data[nan_pix] = 0 - mask[nan_pix] = False - - # If IFU or MOS, mask the fixed-slit area of the image; uses hardwired indexes - if exptype == 'nrs_ifu': - log.info("Masking the fixed slit region for IFU data.") - mask[922:1116, :] = False - elif exptype == 'nrs_msaspec': - # check for any slits defined in the fixed slit quadrant: - # if there is nothing there of interest, mask the whole FS region - slit2msa = input_model.meta.wcs.get_transform('slit_frame', 'msa_frame') - is_fs = [s.quadrant == 5 for s in slit2msa.slits] - if not any(is_fs): - log.info("Masking the fixed slit region for MOS data.") - mask[922:1116, :] = False - else: - log.info("Fixed slits found in MSA definition; " - "not masking the fixed slit region for MOS data.") - - # Mask outliers using sigma clipping stats. - # For BOTS mode, which uses 3D data, loop over each integration separately. - if len(input_model.data.shape) == 3: - for i in range(input_model.data.shape[0]): - _, median, sigma = sigma_clipped_stats(input_model.data[i], mask=~mask[i], mask_value=0, sigma=5.0) - outliers = input_model.data[i] > (median + n_sigma * sigma) - mask[i][outliers] = False - else: - _, median, sigma = sigma_clipped_stats(input_model.data, mask=~mask, mask_value=0, sigma=5.0) - outliers = input_model.data > (median + n_sigma * sigma) - mask[outliers] = False - - - # Return the mask and the record of which pixels were NaN in the input; - # it'll be needed later - return mask, nan_pix - - -def clean_full_frame(detector, image, mask): - """Clean a full-frame (2048x2048) image. - - Parameters - ---------- - detector : str - The name of the detector from which the data originate. - - image : 2D float array - The image to be cleaned. - - mask : 2D bool array - The mask that indicates which pixels are to be used in fitting. - - Returns - ------- - cleaned_image : 2D float array - The cleaned image. - """ - - # Instantiate the cleaner - cleaner = NSClean(detector, mask) - - # Clean the image - try: - cleaned_image = cleaner.clean(image, buff=True) - except np.linalg.LinAlgError: - log.warning("Error cleaning image; step will be skipped") - return None - - return cleaned_image - - -def clean_subarray(detector, image, mask, npix_iter=512, - fc=(1061, 1211, 49943, 49957), - exclude_outliers=True, sigrej=4, minfrac=0.05): - """Clean a subarray image. - - Parameters - ---------- - detector : str - The name of the detector from which the data originate. - - image : 2D float array - The image to be cleaned. - - mask : 2D bool array - The mask that indicates which pixels are to be used in fitting. - - npix_iter : int - Number of pixels to process simultaneously. Default 512. Should - be at least a few hundred to access sub-kHz frequencies in areas - where most pixels are available for fitting. Previous default - behavior corresponds to npix_iter of infinity. - - fc : tuple - Apodizing filter definition. These parameters are tunable. The - defaults happen to work well for NIRSpec BOTS exposures. - 1) Unity gain for f < fc[0] - 2) Cosine roll-off from fc[0] to fc[1] - 3) Zero gain from fc[1] to fc[2] - 4) Cosine roll-on from fc[2] to fc[3] - Default (1061, 1211, 49943, 49957) - - exclude_outliers : bool - Find and mask outliers in the fit? Default True - - sigrej : float - Number of sigma to clip when identifying outliers. Default 4. - - minfrac : float - Minimum fraction of pixels locally available in the mask in - order to attempt a correction. Default 0.05 (i.e., 5%). - - Returns - ------- - cleaned_image : 2D float array - The cleaned image. - """ - - # Flip the image to detector coords. NRS1 requires a transpose - # of the axes, while NRS2 requires a transpose and flip. - if detector == "NRS1": - image = image.transpose() - mask = mask.transpose() - else: - image = image.transpose()[::-1] - mask = mask.transpose()[::-1] - - # We must do the masking of discrepant pixels here: it just - # doesn't work if we wait and do it in the cleaner. This is - # basically copied from lib.py. Use a robust estimator for - # standard deviation, then exclude discrepant pixels and their - # four nearest neighbors from the fit. - - if exclude_outliers: - med = np.median(image[mask]) - std = 1.4825 * np.median(np.abs((image - med)[mask])) - outlier = mask & (np.abs(image - med) > sigrej * std) - - mask = mask&(~outlier) - - # also get four nearest neighbors of flagged pixels - mask[1:] = mask[1:] & (~outlier[:-1]) - mask[:-1] = mask[:-1] & (~outlier[1:]) - mask[:, 1:] = mask[:, 1:] & (~outlier[:, :-1]) - mask[:, :-1] = mask[:, :-1] & (~outlier[:, 1:]) - - # Used to determine the fitting intervals along the slow scan - # direction. Pre-pend a zero so that sum_mask[i] is equal - # to np.sum(mask[:i], axis=1). - - sum_mask = np.array([0] + list(np.cumsum(np.sum(mask, axis=1)))) - - # i1 will be the first row with a nonzero element in the mask - # imax will be the last row with a nonzero element in the mask - - nonzero_mask_element = np.sum(mask, axis=1) > 0 - i1 = np.amin(np.arange(mask.shape[0])[nonzero_mask_element]) - imax = np.amax(np.arange(mask.shape[0])[nonzero_mask_element]) - - i1_vals = [] - di_list = [] - models = [] - while i1 <= imax: - - # Want npix_iter available pixels in this section. If - # there are fewer than 1.5*npix_iter available pixels in - # the rest of the image, just go to the end. - - for k in range(i1 + 1, imax + 2): - if (sum_mask[k] - sum_mask[i1] > npix_iter - and sum_mask[-1] - sum_mask[i1] > 1.5 * npix_iter): - break - - di = k - i1 - - i1_vals += [i1] - di_list += [di] - - # Fit this section only if at least minpct% of the pixels - # are available for finding the background. Don't flag - # outliers section-by-section; we have to do that earlier - # over the full array to get reliable values for the mean - # and standard deviation. - - if np.mean(mask[i1:i1 + di]) > minfrac: - cleaner = NSCleanSubarray(image[i1:i1 + di], mask[i1:i1 + di], - fc=fc, exclude_outliers=False) - models += [cleaner.clean(return_model=True)] - else: - log.warning("Insufficient reference pixels for NSClean around " - "row %d; no correction will be made here." % i1) - models += [np.zeros(image[i1:i1 + di].shape)] - - # If we have reached the end of the array, we are finished. - if k == imax + 1: - break - - # Step forward by half an interval so that we have - # overlapping fitting regions. - - i1 += max(int(np.round(di/2)), 1) - - model = np.zeros(image.shape) - tot_wgt = np.zeros(image.shape) - - # When we combine different corrections computed over - # different intervals, each one the highest weight towards the - # center of its interval and less weight towards the edge. - # Use nonzero weights everywhere so that if only one - # correction is available it gets unit weight when we - # normalize. - - for i in range(len(models)): - wgt = 1.001 - np.abs(np.linspace(-1, 1, di_list[i]))[:, np.newaxis] - model[i1_vals[i]:i1_vals[i] + di_list[i]] += wgt*models[i] - tot_wgt[i1_vals[i]:i1_vals[i] + di_list[i]] += wgt - - # don't divide by zero - tot_wgt[model == 0] = 1 - model /= tot_wgt - cleaned_image = image - model - - # Restore the cleaned image to the science frame - if detector == "NRS1": - cleaned_image = cleaned_image.transpose() - else: - cleaned_image = cleaned_image[::-1].transpose() - - return cleaned_image - - -def do_correction(input_model, mask_spectral_regions, n_sigma, save_mask, user_mask): - - """Apply the NSClean 1/f noise correction - - Parameters - ---------- - input_model : data model object - science data to be corrected - - mask_spectral_regions : bool - Mask slit/slice regions defined in WCS - - n_sigma : float - n-sigma rejection level for finding outliers - - save_mask : bool - switch to indicate whether the mask should be saved - - user_mask : str or None - Path to user-supplied mask image - - Returns - ------- - output_model : `~jwst.datamodel.JwstDataModel` - corrected data - - mask_model : `~jwst.datamodel.JwstDataModel` - pixel mask to be saved or None - """ - - detector = input_model.meta.instrument.detector.upper() - exp_type = input_model.meta.exposure.type - log.info(f'Input exposure type is {exp_type}, detector={detector}') - - output_model = input_model.copy() - - # Check for a user-supplied mask image. If so, use it. - if user_mask is not None: - mask_model = datamodels.open(user_mask) - Mask = (mask_model.data.copy()).astype(np.bool_) - - # Reset and save list of NaN pixels in the input image - nan_pix = np.isnan(input_model.data) - input_model.data[nan_pix] = 0 - Mask[nan_pix] = False - - else: - # Create the pixel mask that'll be used to indicate which pixels - # to include in the 1/f noise measurements. Basically, we're setting - # all illuminated pixels to False, so that they do not get used, and - # setting all unilluminated pixels to True (so they DO get used). - # For BOTS mode the mask will be 3D, to accommodate changes in masked - # pixels per integration. - log.info("Creating mask") - Mask, nan_pix = create_mask(input_model, mask_spectral_regions, n_sigma) - - # Store the mask image in a model, if requested - if save_mask: - if len(Mask.shape) == 3: - mask_model = datamodels.CubeModel(data=Mask) - else: - mask_model = datamodels.ImageModel(data=Mask) - else: - mask_model = None - - log.info(f"Cleaning image {input_model.meta.filename}") - - # Setup for handling 2D or 3D inputs - if len(input_model.data.shape) == 3: - nints = input_model.data.shape[0] - # Check for 3D mask - if len(Mask.shape) == 2: - log.warning("Data are 3D, but mask is 2D. Step will be skipped.") - output_model.meta.cal_step.nsclean = 'SKIPPED' - return output_model, None - else: - nints = 1 - - # Loop over integrations (even if there's only 1) - for i in range(nints): - log.debug(f" working on integration {i+1}") - if len(input_model.data.shape) == 3: - image = np.float32(input_model.data[i]) - mask = Mask[i] - else: - image = np.float32(input_model.data) - mask = Mask - - if input_model.data.shape[-2:] == (2048, 2048): - # Clean a full-frame image - cleaned_image = clean_full_frame(detector, image, mask) - else: - # BOTS and ALLSLITS exposures should be fitting different - # ranges of 1/f frequencies. Be less aggressive with - # fitting higher frequencies in ALLSLITS mode. - if input_model.meta.subarray.name.upper() == "ALLSLITS": - fc = (150, 200, 49943, 49957) - else: - fc = (1061, 1211, 49943, 49957) - - # Clean a subarray image - cleaned_image = clean_subarray(detector, image, mask, fc=fc) - - # Check for failure - if cleaned_image is None: - output_model.meta.cal_step.nsclean = 'SKIPPED' - break - else: - # Store the cleaned image in the output model - if len(output_model.data.shape) == 3: - output_model.data[i] = cleaned_image - else: - output_model.data = cleaned_image - - # Restore NaN's from original image - output_model.data[nan_pix] = np.nan - - # Set completion status - output_model.meta.cal_step.nsclean = 'COMPLETE' - - return output_model, mask_model diff --git a/jwst/nsclean/nsclean_step.py b/jwst/nsclean/nsclean_step.py index 79a61128e1..83c25f42fe 100644 --- a/jwst/nsclean/nsclean_step.py +++ b/jwst/nsclean/nsclean_step.py @@ -1,7 +1,7 @@ from stdatamodels.jwst import datamodels +from jwst.clean_flicker_noise import clean_flicker_noise from ..stpipe import Step -from . import nsclean __all__ = ["NSCleanStep"] @@ -10,15 +10,25 @@ class NSCleanStep(Step): """ NSCleanStep: This step performs 1/f noise correction ("cleaning") of NIRSpec images, using the "NSClean" method. + + NOTE: This step is a deprecated alias to the ``clean_flicker_noise`` step. """ class_alias = "nsclean" spec = """ - mask_spectral_regions = boolean(default=True) # Mask WCS-defined regions + fit_method = option('fft', 'median', default='fft') # Noise fitting algorithm + fit_by_channel = boolean(default=False) # Fit noise separately by amplifier + background_method = option('median', 'model', None, default=None) # Background fitting algorithm + background_box_size = int_list(min=2, max=2, default=None) # Background box size for modeled background + mask_spectral_regions = boolean(default=True) # Mask WCS-defined spectral regions n_sigma = float(default=5.0) # Clipping level for outliers - save_mask = boolean(default=False) # Save the created mask + fit_histogram = boolean(default=False) # Fit a value histogram to derive sigma + single_mask = boolean(default=False) # Make a single mask for all integrations user_mask = string(default=None) # Path to user-supplied mask + save_mask = boolean(default=False) # Save the created mask + save_background = boolean(default=False) # Save the fit background + save_noise = boolean(default=False) # Save the fit noise skip = boolean(default=True) # By default, skip the step """ @@ -29,39 +39,102 @@ def process(self, input): Parameters ---------- input : `~jwst.datamodels.ImageModel`, `~jwst.datamodels.IFUImageModel` - Input datamodel to be corrected + Input datamodel to be corrected. + + fit_method : str, optional + The background fit algorithm to use. Options are 'fft' and 'median'; + 'fft' performs the original NSClean implementation. + + fit_by_channel : bool, optional + If set, flicker noise is fit independently for each detector channel. + Ignored for subarray data and for `fit_method` = 'fft'. + + background_method : {'median', 'model', None} + If 'median', the preliminary background to remove and restore + is a simple median of the background data. If 'model', the + background data is fit with a low-resolution model via + `~photutils.background.Background2D`. If None, the background + value is 0.0. + + background_box_size : tuple of int, optional + Box size for the data grid used by `Background2D` when + `background_method` = 'model'. For best results, use a box size + that evenly divides the input image shape. + + mask_spectral_regions : bool, optional + Mask regions of the image defined by WCS bounding boxes for slits/slices. n_sigma : float, optional - Sigma clipping threshold to be used in detecting outliers in the image + Sigma clipping threshold to be used in detecting outliers in the image. - save_mask : bool, optional - Save the computed mask image + fit_histogram : bool, optional + If set, the 'sigma' used with `n_sigma` for clipping outliers + is derived from a Gaussian fit to a histogram of values. + Otherwise, a simple iterative sigma clipping is performed. + + single_mask : bool, optional + If set, a single mask will be created, regardless of + the number of input integrations. Otherwise, the mask will + be a 3D cube, with one plane for each integration. user_mask : None, str, or `~jwst.datamodels.ImageModel` - Optional user-supplied mask image; path to file or opened datamodel + Optional user-supplied mask image; path to file or opened datamodel. - mask_spectral_regions : bool, optional - Mask regions of the image defined by WCS bounding boxes for slits/slices + save_mask : bool, optional + Save the computed mask image. + + save_background : bool, optional + Save the computed background image. + + save_noise : bool, optional + Save the computed noise image. Returns ------- output_model : `~jwst.datamodels.ImageModel`, `~jwst.datamodels.IFUImageModel` - The 1/f corrected datamodel + The 1/f corrected datamodel. """ + message = ("The 'nsclean' step is a deprecated alias to 'clean_flicker_noise' " + "and will be removed in future builds.") + self.log.warning(message) # Open the input data model with datamodels.open(input) as input_model: # Do the NSClean correction - result = nsclean.do_correction(input_model, self.mask_spectral_regions, self.n_sigma, - self.save_mask, self.user_mask) - output_model, mask_model = result + result = clean_flicker_noise.do_correction( + input_model, self.input_dir, self.fit_method, self.fit_by_channel, + self.background_method, self.background_box_size, + self.mask_spectral_regions, self.n_sigma, self.fit_histogram, + self.single_mask, self.user_mask, + self.save_mask, self.save_background, self.save_noise) + output_model, mask_model, background_model, noise_model, status = result # Save the mask, if requested if self.save_mask and mask_model is not None: - mask_path = self.make_output_path(basepath=input_model.meta.filename, suffix='mask') + mask_path = self.make_output_path( + basepath=input_model.meta.filename, suffix='mask') self.log.info(f"Saving mask file {mask_path}") mask_model.save(mask_path) mask_model.close() + # Save the background, if requested + if self.save_background and background_model is not None: + bg_path = self.make_output_path( + basepath=input_model.meta.filename, suffix='flicker_bkg') + self.log.info(f"Saving background file {bg_path}") + background_model.save(bg_path) + background_model.close() + + # Save the noise, if requested + if self.save_noise and noise_model is not None: + noise_path = self.make_output_path( + basepath=input_model.meta.filename, suffix='flicker_noise') + self.log.info(f"Saving noise file {noise_path}") + noise_model.save(noise_path) + noise_model.close() + + # Set the step completion status + output_model.meta.cal_step.nsclean = status + return output_model diff --git a/jwst/pipeline/calwebb_detector1.py b/jwst/pipeline/calwebb_detector1.py index 16087ee73e..78699ff12f 100644 --- a/jwst/pipeline/calwebb_detector1.py +++ b/jwst/pipeline/calwebb_detector1.py @@ -22,6 +22,7 @@ from ..persistence import persistence_step from ..charge_migration import charge_migration_step from ..jump import jump_step +from ..clean_flicker_noise import clean_flicker_noise_step from ..ramp_fitting import ramp_fit_step from ..gain_scale import gain_scale_step @@ -64,6 +65,7 @@ class Detector1Pipeline(Pipeline): 'persistence': persistence_step.PersistenceStep, 'charge_migration': charge_migration_step.ChargeMigrationStep, 'jump': jump_step.JumpStep, + 'clean_flicker_noise': clean_flicker_noise_step.CleanFlickerNoiseStep, 'ramp_fit': ramp_fit_step.RampFitStep, 'gain_scale': gain_scale_step.GainScaleStep, } @@ -128,6 +130,9 @@ def process(self, input): # apply the jump step input = self.jump(input) + # apply the clean_flicker_noise step + input = self.clean_flicker_noise(input) + # save the corrected ramp data, if requested if self.save_calibrated_ramp: self.save_model(input, 'ramp') diff --git a/jwst/pipeline/calwebb_spec2.py b/jwst/pipeline/calwebb_spec2.py index 37eafd0892..83b542ca51 100644 --- a/jwst/pipeline/calwebb_spec2.py +++ b/jwst/pipeline/calwebb_spec2.py @@ -249,11 +249,11 @@ 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 + # 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'], + calibrated, members_by_type['selfcal'], members_by_type['background'], ) if isinstance(result, datamodels.JwstDataModel): # if step is skipped, unchanged sci exposure is returned diff --git a/jwst/pipeline/clean_noise.cfg b/jwst/pipeline/clean_noise.cfg new file mode 100644 index 0000000000..f528229389 --- /dev/null +++ b/jwst/pipeline/clean_noise.cfg @@ -0,0 +1,2 @@ +name = "clean_flicker_noise" +class = "jwst.clean_flicker_noise.CleanFlickerNoiseStep" diff --git a/jwst/pipeline/tests/test_calwebb_spec2.py b/jwst/pipeline/tests/test_calwebb_spec2.py index d060aa9cd0..71a39232f3 100644 --- a/jwst/pipeline/tests/test_calwebb_spec2.py +++ b/jwst/pipeline/tests/test_calwebb_spec2.py @@ -155,4 +155,4 @@ def test_filenotfounderror_raised(capsys): # Verify the failure is printed to stderr captured = capsys.readouterr() - assert 'FileNotFoundError' in captured.err \ No newline at end of file + assert 'FileNotFoundError' in captured.err diff --git a/jwst/regtest/test_miri_image.py b/jwst/regtest/test_miri_image.py index 5559daf121..192e6dcecb 100644 --- a/jwst/regtest/test_miri_image.py +++ b/jwst/regtest/test_miri_image.py @@ -47,6 +47,26 @@ def run_detector1_with_average_dark_current(rtdata_module): Step.from_cmdline(args) +@pytest.fixture(scope="module") +def run_detector1_with_clean_flicker_noise(rtdata_module): + """Run detector1 pipeline on MIRI imaging data with noise cleaning.""" + rtdata_module.get_data("miri/image/jw01024002001_02101_00001_mirimage_uncal.fits") + + # Run detector1 pipeline only on one of the _uncal files. + # Run optional clean_flicker_noise step, + # saving extra outputs and masking to science regions + args = ["jwst.pipeline.Detector1Pipeline", rtdata_module.input, + "--save_calibrated_ramp=True", + "--steps.clean_flicker_noise.skip=False", + "--steps.clean_flicker_noise.mask_science_regions=True", + "--steps.clean_flicker_noise.save_results=True", + "--steps.clean_flicker_noise.save_mask=True", + "--steps.clean_flicker_noise.save_background=True", + "--steps.clean_flicker_noise.save_noise=True", + ] + Step.from_cmdline(args) + + @pytest.fixture(scope="module") def run_image2(run_detector1, rtdata_module): """Run image2 pipeline on the _rate file, saving intermediate products""" @@ -114,6 +134,29 @@ def test_miri_image_detector1_with_avg_dark_current(run_detector1_with_average_d assert diff.identical, diff.report() +@pytest.mark.bigdata +@pytest.mark.parametrize("suffix", + ["clean_flicker_noise", "mask", + "flicker_bkg", "flicker_noise", + "ramp", "rate", "rateints"]) +def test_miri_image_detector1_with_clean_flicker_noise( + run_detector1_with_clean_flicker_noise, + rtdata_module, fitsdiff_default_kwargs, suffix): + """Test detector1 pipeline for MIRI imaging data with noise cleaning.""" + rtdata = rtdata_module + rtdata.input = "jw01024002001_02101_00001_mirimage_uncal.fits" + output = f"jw01024002001_02101_00001_mirimage_{suffix}.fits" + rtdata.output = output + + rtdata.get_truth(f"truth/test_miri_image_clean_flicker_noise/{output}") + + # Set tolerances so the file comparisons work across architectures + fitsdiff_default_kwargs["rtol"] = 1e-4 + fitsdiff_default_kwargs["atol"] = 1e-4 + diff = FITSDiff(rtdata.output, rtdata.truth, **fitsdiff_default_kwargs) + assert diff.identical, diff.report() + + @pytest.mark.bigdata @pytest.mark.parametrize("suffix", ["assign_wcs", "flat_field", "cal", "i2d"]) def test_miri_image_image2(run_image2, rtdata_module, fitsdiff_default_kwargs, suffix): diff --git a/jwst/regtest/test_nircam_image.py b/jwst/regtest/test_nircam_image.py index 11283baf01..1de07dbf93 100644 --- a/jwst/regtest/test_nircam_image.py +++ b/jwst/regtest/test_nircam_image.py @@ -30,6 +30,24 @@ def run_detector1pipeline(rtdata_module): Step.from_cmdline(args) +@pytest.fixture(scope="module") +def run_detector1_with_clean_flicker_noise(rtdata_module): + """Run detector1 pipeline on NIRCam imaging data with noise cleaning.""" + rtdata_module.get_data("nircam/image/jw01538046001_03105_00001_nrcalong_uncal.fits") + + # Run detector1 pipeline only on one of the _uncal files. + # Run optional clean_flicker_noise step, saving extra outputs + args = ["jwst.pipeline.Detector1Pipeline", rtdata_module.input, + "--save_calibrated_ramp=True", + "--steps.clean_flicker_noise.skip=False", + "--steps.clean_flicker_noise.save_results=True", + "--steps.clean_flicker_noise.save_mask=True", + "--steps.clean_flicker_noise.save_background=True", + "--steps.clean_flicker_noise.save_noise=True", + ] + Step.from_cmdline(args) + + @pytest.fixture(scope="module") def run_image2pipeline(run_detector1pipeline, rtdata_module): """Run calwebb_image2 on NIRCam imaging long data""" @@ -225,3 +243,26 @@ def test_imaging_distortion(rtdata, fitsdiff_default_kwargs): assert_allclose(y, model.meta.wcsinfo.crpix2 - 1) assert_allclose(raout, ra) assert_allclose(decout, dec) + + +@pytest.mark.bigdata +@pytest.mark.parametrize("suffix", + ["clean_flicker_noise", "mask", + "flicker_bkg", "flicker_noise", + "ramp", "rate", "rateints"]) +def test_nircam_image_detector1_with_clean_flicker_noise( + run_detector1_with_clean_flicker_noise, + rtdata_module, fitsdiff_default_kwargs, suffix): + """Test detector1 pipeline for NIRCam imaging data with noise cleaning.""" + rtdata = rtdata_module + rtdata.input = "jw01538046001_03105_00001_nrcalong_uncal.fits" + output = f"jw01538046001_03105_00001_nrcalong_{suffix}.fits" + rtdata.output = output + + rtdata.get_truth(f"truth/test_nircam_image_clean_flicker_noise/{output}") + + # Set tolerances so the file comparisons work across architectures + fitsdiff_default_kwargs["rtol"] = 1e-4 + fitsdiff_default_kwargs["atol"] = 1e-4 + diff = FITSDiff(rtdata.output, rtdata.truth, **fitsdiff_default_kwargs) + assert diff.identical, diff.report() diff --git a/jwst/regtest/test_niriss_image.py b/jwst/regtest/test_niriss_image.py index 7e186986a2..62b268eb6b 100644 --- a/jwst/regtest/test_niriss_image.py +++ b/jwst/regtest/test_niriss_image.py @@ -36,14 +36,46 @@ def run_detector1(rtdata_module): Step.from_cmdline(args) +@pytest.fixture(scope="module") +def run_detector1_with_clean_flicker_noise(rtdata_module): + """Run detector1 pipeline on NIRISS imaging data with noise cleaning.""" + rtdata_module.get_data("niriss/imaging/jw01094001002_02107_00001_nis_uncal.fits") + + # Run detector1 pipeline only on one of the _uncal files. + # Run optional clean_flicker_noise step, + # saving extra outputs and masking to science regions + args = ["jwst.pipeline.Detector1Pipeline", rtdata_module.input, + "--save_calibrated_ramp=True", + "--steps.clean_flicker_noise.skip=False", + "--steps.clean_flicker_noise.save_results=True", + "--steps.clean_flicker_noise.save_mask=True", + "--steps.clean_flicker_noise.save_background=True", + "--steps.clean_flicker_noise.save_noise=True", + ] + Step.from_cmdline(args) + + @pytest.mark.bigdata @pytest.mark.parametrize("suffix", ["dq_init", "saturation", "superbias", "refpix", "linearity", "dark_current", "charge_migration", "jump", "rate", "rateints"]) def test_niriss_image_detector1(run_detector1, rtdata_module, fitsdiff_default_kwargs, suffix): + """Test detector1 pipeline for NIRISS imaging data with noise cleaning. + """ + truth_dir = 'test_niriss_image' + _assert_is_same(rtdata_module, fitsdiff_default_kwargs, suffix, truth_dir) + + +@pytest.mark.bigdata +@pytest.mark.parametrize("suffix", ["clean_flicker_noise", "mask", + "flicker_bkg", "flicker_noise", + "ramp", "rate", "rateints"]) +def test_niriss_image_detector1_with_clean_flicker_noise( + run_detector1_with_clean_flicker_noise, rtdata_module, fitsdiff_default_kwargs, suffix): """Regression test of detector1 pipeline performed on NIRISS imaging data. """ - _assert_is_same(rtdata_module, fitsdiff_default_kwargs, suffix) + truth_dir = 'test_niriss_image_clean_flicker_noise' + _assert_is_same(rtdata_module, fitsdiff_default_kwargs, suffix, truth_dir) @pytest.mark.bigdata @@ -77,14 +109,14 @@ def test_niriss_tweakreg_no_sources(rtdata, fitsdiff_default_kwargs): result.shelve(model, modify=False) -def _assert_is_same(rtdata_module, fitsdiff_default_kwargs, suffix): +def _assert_is_same(rtdata_module, fitsdiff_default_kwargs, suffix, truth_dir): """Assertion helper for the above tests""" rtdata = rtdata_module rtdata.input = "jw01094001002_02107_00001_nis_uncal.fits" output = f"jw01094001002_02107_00001_nis_{suffix}.fits" rtdata.output = output - rtdata.get_truth(f"truth/test_niriss_image/{output}") + rtdata.get_truth(f"truth/{truth_dir}/{output}") # Set tolerances so the crf, rscd and rateints file comparisons work across # architectures diff --git a/jwst/regtest/test_nirspec_irs2_detector1.py b/jwst/regtest/test_nirspec_irs2_detector1.py index 407a8bd2dc..8a8441e26c 100644 --- a/jwst/regtest/test_nirspec_irs2_detector1.py +++ b/jwst/regtest/test_nirspec_irs2_detector1.py @@ -30,6 +30,26 @@ def run_detector1pipeline(rtdata_module): ]) +@pytest.fixture(scope="module") +def run_detector1_with_clean_flicker_noise(rtdata_module): + """Run detector1 pipeline on NIRSpec IRS2 data with noise cleaning.""" + rtdata_module.get_data("nirspec/irs2/jw01335004001_03101_00002_nrs2_uncal.fits") + + # Run detector1 pipeline only on one of the _uncal files. + # Run optional clean_flicker_noise step, + # saving extra outputs and masking to science regions + args = ["jwst.pipeline.Detector1Pipeline", rtdata_module.input, + "--save_calibrated_ramp=True", + "--steps.clean_flicker_noise.skip=False", + "--steps.clean_flicker_noise.mask_science_regions=True", + "--steps.clean_flicker_noise.save_results=True", + "--steps.clean_flicker_noise.save_mask=True", + "--steps.clean_flicker_noise.save_background=True", + "--steps.clean_flicker_noise.save_noise=True", + ] + Step.from_cmdline(args) + + @pytest.mark.bigdata @pytest.mark.parametrize("suffix", ['dq_init', 'saturation', 'superbias', 'refpix', 'linearity', 'dark_current', 'jump', @@ -51,3 +71,24 @@ def test_nirspec_irs2_detector1(run_detector1pipeline, rtdata_module, fitsdiff_default_kwargs["atol"] = 1e-4 diff = FITSDiff(rtdata.output, rtdata.truth, **fitsdiff_default_kwargs) assert diff.identical, diff.report() + + +@pytest.mark.bigdata +@pytest.mark.parametrize("suffix", ["clean_flicker_noise", "mask", + "flicker_bkg", "flicker_noise", + "ramp", "rate", "rateints"]) +def test_nirspec_irs2_detector1_with_clean_flicker_noise( + run_detector1_with_clean_flicker_noise, rtdata_module, + fitsdiff_default_kwargs, suffix): + """Test detector1 pipeline for NIRSpec IRS2 data with noise cleaning.""" + rtdata = rtdata_module + + output_filename = f"jw01335004001_03101_00002_nrs2_{suffix}.fits" + rtdata.output = output_filename + rtdata.get_truth(f"truth/test_nirspec_irs2_clean_flicker_noise/{output_filename}") + + # Set tolerances so comparisons work across architectures + fitsdiff_default_kwargs["rtol"] = 1e-4 + fitsdiff_default_kwargs["atol"] = 1e-4 + diff = FITSDiff(rtdata.output, rtdata.truth, **fitsdiff_default_kwargs) + assert diff.identical, diff.report() diff --git a/jwst/step.py b/jwst/step.py index 1c9b67b58f..fe051a7b62 100644 --- a/jwst/step.py +++ b/jwst/step.py @@ -7,6 +7,7 @@ from .badpix_selfcal.badpix_selfcal_step import BadpixSelfcalStep from .barshadow.barshadow_step import BarShadowStep from .charge_migration.charge_migration_step import ChargeMigrationStep +from .clean_flicker_noise.clean_flicker_noise_step import CleanFlickerNoiseStep from .combine_1d.combine_1d_step import Combine1dStep from .coron.stack_refs_step import StackRefsStep from .coron.align_refs_step import AlignRefsStep @@ -76,6 +77,7 @@ "AlignRefsStep", "KlipStep", "HlspStep", + "CleanFlickerNoiseStep", "CubeBuildStep", "CubeSkyMatchStep", "DarkCurrentStep", diff --git a/jwst/stpipe/integration.py b/jwst/stpipe/integration.py index 597a6afe3e..7b787b47a7 100644 --- a/jwst/stpipe/integration.py +++ b/jwst/stpipe/integration.py @@ -44,6 +44,7 @@ def get_steps(): ("jwst.step.AlignRefsStep", 'align_refs', False), ("jwst.step.KlipStep", 'klip', False), ("jwst.step.HlspStep", 'hlsp', False), + ("jwst.step.CleanFlickerNoiseStep", 'clean_flicker_noise', False), ("jwst.step.CubeBuildStep", 'cube_build', False), ("jwst.step.CubeSkyMatchStep", 'cube_skymatch', False), ("jwst.step.DarkCurrentStep", 'dark_current', False),