diff --git a/CHANGES.rst b/CHANGES.rst index b6add1e8a7..1346d55e59 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -46,7 +46,9 @@ associations - Added default values for new non-header keywords (``MOSTILNO`` and ``DITHPTIN``) to the default values in the ``asn_make_pool`` script. [#8508] -- Create WFSS Pure-Parallel associations [#8528] +- Create WFSS Pure-Parallel associations. [#8528] + +- Add NIRSpec optical path constraints for TSO associations. [#8537] combine_1d ---------- @@ -245,6 +247,9 @@ pixel_replace - Moved pixel_replace in the calwebb_spec2 pipeline and added it to the calwebb_spec3 pipeline. In both pipelines it is now executed immediately before resample_spec/cube_build. [#8409] +- Added estimated errors and variances for replaced pixels, following the + interpolation scheme used for the data. [#8504] + ramp_fitting ------------ @@ -274,6 +279,8 @@ resample - Change `fillval` parameter default from INDEF to NaN [#8488] +- Removed the use of the `drizpars` reference file [#8546] + resample_spec ------------- @@ -285,6 +292,8 @@ resample_spec when the slit was closely aligned with the RA direction sky. [#8511] +- Removed the use of the `drizpars` reference file [#8546] + residual_fringe --------------- diff --git a/docs/jwst/pixel_replace/arguments.rst b/docs/jwst/pixel_replace/arguments.rst index aaaf2a7cd3..aeebc1711d 100644 --- a/docs/jwst/pixel_replace/arguments.rst +++ b/docs/jwst/pixel_replace/arguments.rst @@ -5,11 +5,12 @@ The ``pixel_replace`` step has the following step-specific arguments: ``--algorithm`` (str, default='fit_profile') This sets the method used to estimate flux values for bad pixels. The default 'fit_profile' uses a profile - fit to adjacent column values. The minimum gradient ('mingrad') method is also available for the MIRI MRS. + fit to adjacent column values. The minimum gradient ('mingrad') method, developed for MIRI MRS data, may + also be used for all instruments and modes. ``--n_adjacent_cols`` (int, default=3) Number of adjacent columns (on either side of column containing a bad pixel) to use in creation of the source profile, in cross-dispersion direction. The total number of columns used in the profile will be twice this number; on array edges, the total number of columns contributing to the source profile will be less than ``2 * n_adjacent_cols``. - + Ignored when ``algorithm = 'mingrad'``. diff --git a/docs/jwst/pixel_replace/main.rst b/docs/jwst/pixel_replace/main.rst index 8c2c34078b..e939c5a218 100644 --- a/docs/jwst/pixel_replace/main.rst +++ b/docs/jwst/pixel_replace/main.rst @@ -15,6 +15,8 @@ flagged as ``DO_NOT_USE`` in 2-D extracted spectra using interpolation methods, prior to rectification in the :ref:`resample_spec ` step. ``pixel_replace`` inserts these estimates into the 2-D data array, unsets the ``DO_NOT_USE`` flag, and sets the ``FLUX_ESTIMATED`` flag for each affected pixel. +Error values and variance components for the replaced pixels are similarly updated with +estimated values, following the same interpolation method as is used for the data. This step is provided as a cosmetic feature and, for that reason, should be used with caution. diff --git a/docs/jwst/pixel_replace/reference_files.rst b/docs/jwst/pixel_replace/reference_files.rst index 6e51bd5833..69605419f7 100644 --- a/docs/jwst/pixel_replace/reference_files.rst +++ b/docs/jwst/pixel_replace/reference_files.rst @@ -1,4 +1,4 @@ Reference File ============== -This step does not use any reference file. +This step does not use any reference files. diff --git a/docs/jwst/references_general/drizpars_reffile.inc b/docs/jwst/references_general/drizpars_reffile.inc deleted file mode 100644 index faeb08009f..0000000000 --- a/docs/jwst/references_general/drizpars_reffile.inc +++ /dev/null @@ -1,51 +0,0 @@ -.. _drizpars_reffile: - -DRIZPARS Reference File ------------------------ - -:REFTYPE: DRIZPARS -:Data model: `~jwst.datamodels.DrizParsModel` - -The DRIZPARS reference file contains various drizzle parameter values that -control the characteristics of a drizzled image and how it is built. - -.. include:: ../references_general/drizpars_selection.inc - -.. include:: ../includes/standard_keywords.inc - -Type Specific Keywords for DRIZPARS -+++++++++++++++++++++++++++++++++++++ -No additional specific keywords are required in DRIZPARS reference -files, because CRDS selection is based only on the instrument name -(see :ref:`drizpars_selectors`). - -Reference File Format -+++++++++++++++++++++ -DRIZPARS reference files are FITS format, with 1 BINTABLE extension. -The FITS primary HDU does not contain a data array. -The format and content of the file is as follows: - -======== ======== ===== ======================= ========= -EXTNAME XTENSION NAXIS Dimensions Data type -======== ======== ===== ======================= ========= -DRIZPARS BINTABLE 2 TFIELDS = 7 N/A -======== ======== ===== ======================= ========= - -The DRIZPARS extension contains various step parameter values to be -used when processing certain types of image collections. The first -two columns (numimages and filter) are used as row selectors within -the table. Image collections that match those selectors then use the -parameter values specified in the remainder of that table row. -The table contains the following 7 columns: - -=========== ======= =========================================== -TTYPE TFORM Description -=========== ======= =========================================== -numimages integer The number of images to be combined -filter integer The filter used to obtain the images -pixfrac float The pixel "shrinkage" fraction -kernel string The kernel function used to distribute flux -fillval float Value assigned to pixels with no input flux -wht_type string The input image weighting type -stepsize integer Output WCS grid interpolation step size -=========== ======= =========================================== diff --git a/docs/jwst/references_general/drizpars_selection.inc b/docs/jwst/references_general/drizpars_selection.inc deleted file mode 100644 index d45fd71e64..0000000000 --- a/docs/jwst/references_general/drizpars_selection.inc +++ /dev/null @@ -1,16 +0,0 @@ -.. _drizpars_selectors: - -Reference Selection Keywords for DRIZPARS -+++++++++++++++++++++++++++++++++++++++++ -CRDS selects appropriate DRIZPARS references based on the following keywords. -DRIZPARS is not applicable for instruments not in the table. -All keywords used for file selection are *required*. - -========== ======== -Instrument Keywords -========== ======== -MIRI INSTRUME -NIRCam INSTRUME -NIRISS INSTRUME -========== ======== - diff --git a/docs/jwst/references_general/references_general.rst b/docs/jwst/references_general/references_general.rst index 417b6fd894..550eb4dad8 100644 --- a/docs/jwst/references_general/references_general.rst +++ b/docs/jwst/references_general/references_general.rst @@ -150,8 +150,6 @@ documentation on each reference file. +-----------------------------------------------+--------------------------------------------------+ | :ref:`refpix ` | :ref:`REFPIX ` | +-----------------------------------------------+--------------------------------------------------+ -| :ref:`resample ` | :ref:`DRIZPARS ` | -+-----------------------------------------------+--------------------------------------------------+ | :ref:`reset ` | :ref:`RESET ` | +-----------------------------------------------+--------------------------------------------------+ | :ref:`residual_fringe ` | :ref:`FRINGEFREQ ` | @@ -204,8 +202,6 @@ documentation on each reference file. +--------------------------------------------------+-----------------------------------------------+ | :ref:`DISTORTION ` | :ref:`assign_wcs ` | +--------------------------------------------------+-----------------------------------------------+ -| :ref:`DRIZPARS ` | :ref:`resample ` | -+--------------------------------------------------+-----------------------------------------------+ | :ref:`EMICORR ` | :ref:`emicorr ` | +--------------------------------------------------+-----------------------------------------------+ | :ref:`EXTRACT1D ` | :ref:`extract_1d ` | diff --git a/docs/jwst/resample/main.rst b/docs/jwst/resample/main.rst index 08e0069f2d..9291519018 100644 --- a/docs/jwst/resample/main.rst +++ b/docs/jwst/resample/main.rst @@ -15,11 +15,9 @@ The ``resample`` step can take as input either: #. a single 2D input image #. an association table (in json format) -The defined parameters for the drizzle operation itself get -provided by the DRIZPARS reference file (from CRDS). The exact values -used depends on the number of input images being combined and the filter -being used. Other information may be added as selection criteria later, -but for now, only basic information is used. +The parameters for the drizzle operation are provided via ``resample`` +step parameters, which can be overridden by a step parameter reference +file from CRDS. The output product gets defined using the WCS information of all inputs, even if it is just a single input image. The output WCS defines a diff --git a/docs/jwst/resample/reference_files.rst b/docs/jwst/resample/reference_files.rst index d292c96f09..75ad1e0057 100644 --- a/docs/jwst/resample/reference_files.rst +++ b/docs/jwst/resample/reference_files.rst @@ -1,5 +1,3 @@ Reference File ============== -The ``resample`` step uses the DRIZPARS reference file. - -.. include:: ../references_general/drizpars_reffile.inc +The ``resample`` step does not use any reference files. diff --git a/jwst/associations/lib/rules_level2b.py b/jwst/associations/lib/rules_level2b.py index f176a3ffc4..0ddb850d23 100644 --- a/jwst/associations/lib/rules_level2b.py +++ b/jwst/associations/lib/rules_level2b.py @@ -474,6 +474,24 @@ def __init__(self, *args, **kwargs): ) ], reduce=Constraint.notany + ), + # Don't allow NIRSpec invalid optical paths in spec2 + Constraint( + [ + Constraint([ + DMSAttrConstraint( + name='exp_type', + sources=['exp_type'], + value=('nrs_brightobj') + ), + SimpleConstraint( + value=False, + test=lambda value, item: nrsfss_valid_detector(item) == value, + force_unique=False + ), + ]), + ], + reduce=Constraint.notany ) ]) diff --git a/jwst/associations/lib/rules_level3.py b/jwst/associations/lib/rules_level3.py index 1e93204376..ca6db69a5e 100644 --- a/jwst/associations/lib/rules_level3.py +++ b/jwst/associations/lib/rules_level3.py @@ -881,19 +881,19 @@ def __init__(self, *args, **kwargs): DMSAttrConstraint( name='restricted_grism', sources=['exp_type'], - value = ('nrc_tsgrism') + value=('nrc_tsgrism') ), DMSAttrConstraint( name='grism_clear', sources=['pupil'], - value = ('clear') + value=('clear') ), ]), Constraint([ DMSAttrConstraint( name='restricted_ts', sources=['exp_type'], - value = 'nrc_tsgrism' + value='nrc_tsgrism' ), DMSAttrConstraint( name='module', @@ -911,12 +911,30 @@ def __init__(self, *args, **kwargs): DMSAttrConstraint( name='exp_type', sources=['exp_type'], - value = ('nis_soss') + value=('nis_soss') ), DMSAttrConstraint( name='nints', sources=['nints'], - value = ('1') + value=('1') + ), + ]), + ], + reduce=Constraint.notany + ), + # Don't allow NIRSpec invalid optical paths in TSO3 + Constraint( + [ + Constraint([ + DMSAttrConstraint( + name='exp_type', + sources=['exp_type'], + value=('nrs_brightobj') + ), + SimpleConstraint( + value=False, + test=lambda value, item: nrsfss_valid_detector(item) == value, + force_unique=False ), ]), ], diff --git a/jwst/associations/tests/data/pool_021_tso.csv b/jwst/associations/tests/data/pool_021_tso.csv index 0850c339c9..bb413cbce2 100644 --- a/jwst/associations/tests/data/pool_021_tso.csv +++ b/jwst/associations/tests/data/pool_021_tso.csv @@ -17,7 +17,7 @@ set acid|OBS_ID|PROGRAM|OBS_NUM|VISIT|VISIT_ID|VISITGRP|WFSVISIT|SEQ_ID|ACT_ID|E @!fmt_fname(next(expnum))|V99009003001P0000000002102|99009|@!next(obsnum)|1|99009003001|2|NULL|1|2|1|NIS_SOSS|NULL|NULL|NIRISS|CENTRAL|NULL|1|M-33|NULL|FIXED|NIRISS WFSS|SCIENCE|2|1|1|1|1|NULL|F140M|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NIS_CEN|2015_1|FULL|NULL|NULL|NULL|@!fmt_cand([(next(obsnum), 'OBSERVATION')])|N|N|N|T @!fmt_fname(next(expnum))|V99009003001P0000000002102|99009|@!next(obsnum)|1|99009003001|2|NULL|1|2|1|NRC_TSIMAGE|NULL|NULL|NIRCAM|CENTRAL|NULL|1|M-33|NULL|FIXED|NIRISS WFSS|SCIENCE|2|1|1|1|1|NULL|F140M|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NIS_CEN|2015_1|FULL|NULL|NULL|NULL|@!fmt_cand([(next(obsnum), 'OBSERVATION')])|N|N|N|T @!fmt_fname(next(expnum))|V99009003001P0000000002102|99009|@!next(obsnum)|1|99009003001|2|NULL|1|2|1|NRC_TSGRISM|NULL|NULL|NIRCAM|CENTRAL|NULL|1|M-33|NULL|FIXED|NIRISS WFSS|SCIENCE|2|1|1|1|1|NULL|F140M|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NIS_CEN|2015_1|FULL|NULL|NULL|NULL|@!fmt_cand([(next(obsnum), 'OBSERVATION')])|N|N|N|T -@!fmt_fname(next(expnum))|V99009003001P0000000002102|99009|@!next(obsnum)|1|99009003001|2|NULL|1|2|1|NRS_BRIGHTOBJ|NULL|NULL|NIRSPEC|CENTRAL|NULL|1|M-33|NULL|FIXED|NIRISS WFSS|SCIENCE|2|1|1|1|1|NULL|F140M|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NIS_CEN|2015_1|FULL|NULL|NULL|NULL|@!fmt_cand([(next(obsnum), 'OBSERVATION')])|N|N|N|T +@!fmt_fname(next(expnum))|V99009003001P0000000002102|99009|@!next(obsnum)|1|99009003001|2|NULL|1|2|1|NRS_BRIGHTOBJ|NULL|NULL|NIRSPEC|NRS1|NULL|1|M-33|NULL|FIXED|NIRISS WFSS|SCIENCE|2|1|1|1|1|NULL|CLEAR|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NIS_CEN|2015_1|FULL|PRISM|S1600A1|NULL|@!fmt_cand([(next(obsnum), 'OBSERVATION')])|N|N|N|T # # Set ACID to how many have been used. set acid|OBS_ID|PROGRAM|OBS_NUM|VISIT|VISIT_ID|VISITGRP|WFSVISIT|SEQ_ID|ACT_ID|EXPOSURE|EXP_TYPE|NEXPOSUR|EXPCOUNT|INSTRUME|DETECTOR|CHANNEL|TARGETID|TARGPROP|TARGNAME|TARGTYPE|TEMPLATE|PNTGTYPE|PNTG_SEQ|TARGORDN|EXPSPCIN|DITHPTIN|MOSTILNO|MODULE|FILTER|PUPIL|DITHERID|PATTTYPE|PATTSTRT|NUMDTHPT|PATTSIZE|SUBPXPNS|PATT_NUM|SUBPXNUM|SUBPIXEL|APERNAME|SDP_VER|SUBARRAY|GRATING|FXD_SLIT|BAND|@!acid.set(acid.value)|BACKGROUND|IS_IMPRT|IS_PSF|IS_TSO diff --git a/jwst/associations/tests/data/pool_022_tso_noflag.csv b/jwst/associations/tests/data/pool_022_tso_noflag.csv index d49d33ec36..afa7a192ab 100644 --- a/jwst/associations/tests/data/pool_022_tso_noflag.csv +++ b/jwst/associations/tests/data/pool_022_tso_noflag.csv @@ -17,7 +17,7 @@ set acid|OBS_ID|PROGRAM|OBS_NUM|VISIT|VISIT_ID|VISITGRP|WFSVISIT|SEQ_ID|ACT_ID|E @!fmt_fname(next(expnum))|V99009003001P0000000002102|99009|@!next(obsnum)|1|99009003001|2|NULL|1|2|1|NIS_SOSS|NULL|NULL|NIRISS|CENTRAL|NULL|1|M-33|NULL|FIXED|NIRISS WFSS|SCIENCE|2|1|1|1|1|NULL|F140M|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NIS_CEN|2015_1|FULL|NULL|NULL|NULL|@!fmt_cand([(next(obsnum), 'OBSERVATION')])|N|N|N|T @!fmt_fname(next(expnum))|V99009003001P0000000002102|99009|@!next(obsnum)|1|99009003001|2|NULL|1|2|1|NRC_TSIMAGE|NULL|NULL|NIRCAM|CENTRAL|NULL|1|M-33|NULL|FIXED|NIRISS WFSS|SCIENCE|2|1|1|1|1|NULL|F140M|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NIS_CEN|2015_1|FULL|NULL|NULL|NULL|@!fmt_cand([(next(obsnum), 'OBSERVATION')])|N|N|N|T @!fmt_fname(next(expnum))|V99009003001P0000000002102|99009|@!next(obsnum)|1|99009003001|2|NULL|1|2|1|NRC_TSGRISM|NULL|NULL|NIRCAM|CENTRAL|NULL|1|M-33|NULL|FIXED|NIRISS WFSS|SCIENCE|2|1|1|1|1|NULL|F140M|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NIS_CEN|2015_1|FULL|NULL|NULL|NULL|@!fmt_cand([(next(obsnum), 'OBSERVATION')])|N|N|N|T -@!fmt_fname(next(expnum))|V99009003001P0000000002102|99009|@!next(obsnum)|1|99009003001|2|NULL|1|2|1|NRS_BRIGHTOBJ|NULL|NULL|NIRSPEC|CENTRAL|NULL|1|M-33|NULL|FIXED|NIRISS WFSS|SCIENCE|2|1|1|1|1|NULL|F140M|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NIS_CEN|2015_1|FULL|NULL|NULL|NULL|@!fmt_cand([(next(obsnum), 'OBSERVATION')])|N|N|N|T +@!fmt_fname(next(expnum))|V99009003001P0000000002102|99009|@!next(obsnum)|1|99009003001|2|NULL|1|2|1|NRS_BRIGHTOBJ|NULL|NULL|NIRSPEC|NRS1|NULL|1|M-33|NULL|FIXED|NIRISS WFSS|SCIENCE|2|1|1|1|1|NULL|CLEAR|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NULL|NIS_CEN|2015_1|FULL|PRISM|S1600A1|NULL|@!fmt_cand([(next(obsnum), 'OBSERVATION')])|N|N|N|T # # Set ACID to how many have been used. set acid|OBS_ID|PROGRAM|OBS_NUM|VISIT|VISIT_ID|VISITGRP|WFSVISIT|SEQ_ID|ACT_ID|EXPOSURE|EXP_TYPE|NEXPOSUR|EXPCOUNT|INSTRUME|DETECTOR|CHANNEL|TARGETID|TARGPROP|TARGNAME|TARGTYPE|TEMPLATE|PNTGTYPE|PNTG_SEQ|TARGORDN|EXPSPCIN|DITHPTIN|MOSTILNO|MODULE|FILTER|PUPIL|DITHERID|PATTTYPE|PATTSTRT|NUMDTHPT|PATTSIZE|SUBPXPNS|PATT_NUM|SUBPXNUM|SUBPIXEL|APERNAME|SDP_VER|SUBARRAY|GRATING|FXD_SLIT|BAND|@!acid.set(acid.value)|BACKGROUND|IS_IMPRT|IS_PSF|TSOVISIT diff --git a/jwst/msaflagopen/msaflagopen_step.py b/jwst/msaflagopen/msaflagopen_step.py index afe5f7b1b7..9e2cf3ffc1 100755 --- a/jwst/msaflagopen/msaflagopen_step.py +++ b/jwst/msaflagopen/msaflagopen_step.py @@ -1,5 +1,7 @@ from stdatamodels.jwst import datamodels +from stpipe.crds_client import reference_uri_to_cache_path + from ..stpipe import Step from . import msaflag_open @@ -51,7 +53,6 @@ def process(self, input): def create_reference_filename_dictionary(input_model): reffiles = {} - a = Step() reffiles['distortion'] = input_model.meta.ref_file.distortion.name reffiles['filteroffset'] = input_model.meta.ref_file.filteroffset.name reffiles['specwcs'] = input_model.meta.ref_file.filteroffset.name @@ -72,5 +73,5 @@ def create_reference_filename_dictionary(input_model): # Convert from crds protocol to absolute filenames for key in reffiles.keys(): if reffiles[key].startswith('crds://'): - reffiles[key] = a.reference_uri_to_cache_path(reffiles[key], input_model.crds_observatory) + reffiles[key] = reference_uri_to_cache_path(reffiles[key], input_model.crds_observatory) return reffiles diff --git a/jwst/pipeline/calwebb_spec3.py b/jwst/pipeline/calwebb_spec3.py index 1eeafffd62..9357903a17 100644 --- a/jwst/pipeline/calwebb_spec3.py +++ b/jwst/pipeline/calwebb_spec3.py @@ -227,10 +227,11 @@ def process(self, input): for cal_array in result: cal_array.meta.asn.table_name = op.basename(input_models.asn_table_name) result = self.outlier_detection(result) + # interpolate pixels that have a NaN value or are flagged # as DO_NOT_USE or NON_SCIENCE. - print(result) result = self.pixel_replace(result) + # Resample time. Dependent on whether the data is IFU or not. resample_complete = None if exptype in IFU_EXPTYPES: diff --git a/jwst/pixel_replace/pixel_replace.py b/jwst/pixel_replace/pixel_replace.py index bdd96f2bbd..75d78131aa 100644 --- a/jwst/pixel_replace/pixel_replace.py +++ b/jwst/pixel_replace/pixel_replace.py @@ -71,6 +71,7 @@ def replace(self): to each 2D spectrum in input. """ # ImageModel inputs (MIR_LRS-FIXEDSLIT) + # or 2D SlitModel inputs (e.g. NRS_FIXEDSLIT in spec3) if (isinstance(self.input, datamodels.ImageModel) or (isinstance(self.input, datamodels.SlitModel) and self.input.data.ndim == 2)): @@ -78,94 +79,105 @@ def replace(self): n_replaced = np.count_nonzero(self.output.dq & self.FLUX_ESTIMATED) log.info(f"Input model had {n_replaced} pixels replaced.") elif isinstance(self.input, datamodels.IFUImageModel): - # Attempt to run pixel replacement on each throw of the IFU slicer - # individually. - xx, yy = np.indices(self.input.data.shape) - if self.input.meta.exposure.type == 'MIR_MRS': - # Mingrad method - if (self.pars['algorithm'] == 'mingrad'): - new_model = self.algorithm(self.input) - self.output = new_model + # Attempt to run pixel replacement on each throw of the IFU slicer + # individually. + xx, yy = np.indices(self.input.data.shape) + if self.input.meta.exposure.type == 'MIR_MRS': + if self.pars['algorithm'] == 'mingrad': + # mingrad method + new_model = self.algorithm(self.input) + self.output = new_model + else: # fit_profile method - else: - _, beta_array, _ = self.input.meta.wcs.transform('detector', 'alpha_beta', yy, xx) - unique_beta = np.unique(beta_array) - unique_beta = unique_beta[~np.isnan(unique_beta)] - for i, beta in enumerate(unique_beta): - # Define a mask that is True where this trace is located - trace_mask = (beta_array == beta) - trace_model = self.input.copy() - trace_model.dq = np.where( - # When not in this trace, set NON_SCIENCE and DO_NOT_USE - ~trace_mask, - trace_model.dq | self.DO_NOT_USE | self.NON_SCIENCE, - trace_model.dq - ) - - trace_model = self.algorithm(trace_model) - self.output.data = np.where( - # Where trace is located, set replaced values - trace_mask, - trace_model.data, - self.output.data - ) - - self.output.dq = np.where( - # Where trace is located, set replaced values - trace_mask, - trace_model.dq, - self.output.dq - ) - - n_replaced = np.count_nonzero(trace_model.dq & self.FLUX_ESTIMATED) - log.info(f"Input MRS frame had {n_replaced} pixels replaced in IFU slice {i+1}.") - trace_model.close() - - n_replaced = np.count_nonzero(self.output.dq & self.FLUX_ESTIMATED) - log.info(f"Input MRS frame had {n_replaced} total pixels replaced.") + _, beta_array, _ = self.input.meta.wcs.transform('detector', 'alpha_beta', yy, xx) + unique_beta = np.unique(beta_array) + unique_beta = unique_beta[~np.isnan(unique_beta)] + for i, beta in enumerate(unique_beta): + # Define a mask that is True where this trace is located + trace_mask = (beta_array == beta) + trace_model = self.input.copy() + trace_model.dq = np.where( + # When not in this trace, set NON_SCIENCE and DO_NOT_USE + ~trace_mask, + trace_model.dq | self.DO_NOT_USE | self.NON_SCIENCE, + trace_model.dq + ) + + trace_model = self.algorithm(trace_model) + self.output.data = np.where( + # Where trace is located, set replaced values + trace_mask, + trace_model.data, + self.output.data + ) + + # do the same for dq, err, and var + self.output.dq = np.where( + trace_mask, trace_model.dq, self.output.dq) + self.output.err = np.where( + trace_mask, trace_model.err, self.output.err) + self.output.var_poisson = np.where( + trace_mask, trace_model.var_poisson, self.output.var_poisson) + self.output.var_rnoise = np.where( + trace_mask, trace_model.var_rnoise, self.output.var_rnoise) + self.output.var_flat = np.where( + trace_mask, trace_model.var_flat, self.output.var_flat) + + n_replaced = np.count_nonzero(trace_model.dq & self.FLUX_ESTIMATED) + log.info(f"Input MRS frame had {n_replaced} pixels replaced in IFU slice {i+1}.") + + trace_model.close() + + n_replaced = np.count_nonzero(self.output.dq & self.FLUX_ESTIMATED) + log.info(f"Input MRS frame had {n_replaced} total pixels replaced.") + else: + if self.pars['algorithm'] == 'mingrad': + # mingrad method + new_model = self.algorithm(self.input) + self.output = new_model else: - # Mingrad method - if (self.pars['algorithm'] == 'mingrad'): - new_model = self.algorithm(self.input) - self.output = new_model # fit_profile method - iterate over IFU slices - else: - for i in range(30): - slice_wcs = nrs_wcs_set_input(self.input, i) - _, _, wave = slice_wcs.transform('detector', 'slicer', yy, xx) - - # Define a mask that is True where this trace is located - trace_mask = (wave > 0) - trace_model = self.input.copy() - trace_model.dq = np.where( - # When not in this trace, set NON_SCIENCE and DO_NOT_USE - ~trace_mask, - trace_model.dq | self.DO_NOT_USE | self.NON_SCIENCE, - trace_model.dq - ) - - trace_model = self.algorithm(trace_model) - self.output.data = np.where( - # Where trace is located, set replaced values - trace_mask, - trace_model.data, - self.output.data - ) - - self.output.dq = np.where( - # Where trace is located, set replaced values - trace_mask, - trace_model.dq, - self.output.dq - ) - - n_replaced = np.count_nonzero(trace_model.dq & self.FLUX_ESTIMATED) - log.info(f"Input NRS_IFU frame had {n_replaced} pixels replaced in IFU slice {i + 1}.") - - trace_model.close() - - n_replaced = np.count_nonzero(self.output.dq & self.FLUX_ESTIMATED) - log.info(f"Input NRS_IFU frame had {n_replaced} total pixels replaced.") + for i in range(30): + slice_wcs = nrs_wcs_set_input(self.input, i) + _, _, wave = slice_wcs.transform('detector', 'slicer', yy, xx) + + # Define a mask that is True where this trace is located + trace_mask = (wave > 0) + trace_model = self.input.copy() + trace_model.dq = np.where( + # When not in this trace, set NON_SCIENCE and DO_NOT_USE + ~trace_mask, + trace_model.dq | self.DO_NOT_USE | self.NON_SCIENCE, + trace_model.dq + ) + + trace_model = self.algorithm(trace_model) + self.output.data = np.where( + # Where trace is located, set replaced values + trace_mask, + trace_model.data, + self.output.data + ) + + # do the same for dq, err, and var + self.output.dq = np.where( + trace_mask, trace_model.dq, self.output.dq) + self.output.err = np.where( + trace_mask, trace_model.err, self.output.err) + self.output.var_poisson = np.where( + trace_mask, trace_model.var_poisson, self.output.var_poisson) + self.output.var_rnoise = np.where( + trace_mask, trace_model.var_rnoise, self.output.var_rnoise) + self.output.var_flat = np.where( + trace_mask, trace_model.var_flat, self.output.var_flat) + + n_replaced = np.count_nonzero(trace_model.dq & self.FLUX_ESTIMATED) + log.info(f"Input NRS_IFU frame had {n_replaced} pixels replaced in IFU slice {i + 1}.") + + trace_model.close() + + n_replaced = np.count_nonzero(self.output.dq & self.FLUX_ESTIMATED) + log.info(f"Input NRS_IFU frame had {n_replaced} total pixels replaced.") # MultiSlitModel inputs (WFSS, NRS_FIXEDSLIT, ?) elif isinstance(self.input, datamodels.MultiSlitModel): @@ -180,21 +192,30 @@ def replace(self): self.output.slits[i] = slit_replaced - # CubeModel inputs are TSO (so far?); only SlitModel seen so far is NRS_BRIGHTOBJ, also requiring - # a re-packaging of the data into 2D inputs for the algorithm. + # CubeModel inputs are TSO (so far?); SlitModel may be NRS_BRIGHTOBJ, + # also requiring a re-packaging of the data into 2D inputs for the algorithm elif isinstance(self.input, (datamodels.CubeModel, datamodels.SlitModel)): - # Initial attempt looped over model.meta.exposure.nints, but test data had mismatch. Could change this. for i in range(len(self.input.data)): - dummy_model = datamodels.ImageModel(data=self.input.data[i], dq=self.input.dq[i]) - dummy_model.update(self.input) - dummy_replaced = self.algorithm(dummy_model) - n_replaced = np.count_nonzero(dummy_replaced.dq & self.FLUX_ESTIMATED) + img_model = datamodels.ImageModel( + data=self.input.data[i], dq=self.input.dq[i], + err=self.input.err[i], + var_poisson=self.input.var_poisson[i], + var_rnoise=self.input.var_rnoise[i], + var_flat=self.input.var_flat[i], + ) + img_model.update(self.input) + img_replaced = self.algorithm(img_model) + n_replaced = np.count_nonzero(img_replaced.dq & self.FLUX_ESTIMATED) log.info(f"Input TSO integration {i} had {n_replaced} pixels replaced.") - self.output.data[i] = dummy_replaced.data - self.output.dq[i] = dummy_replaced.dq - dummy_replaced.close() - dummy_model.close() + self.output.data[i] = img_replaced.data + self.output.dq[i] = img_replaced.dq + self.output.err[i] = img_replaced.err + self.output.var_poisson[i] = img_replaced.var_poisson + self.output.var_rnoise[i] = img_replaced.var_rnoise + self.output.var_flat[i] = img_replaced.var_flat + img_replaced.close() + img_model.close() else: # This should never happen, as these should be caught in the step code. @@ -203,10 +224,16 @@ def replace(self): def fit_profile(self, model): """ + Replace pixels with the profile fit method. + Fit a profile to adjacent columns, scale profile to column with missing pixel(s), and find flux estimate from scaled profile. + Error and variance values for the replaced pixels + are similarly estimated, using the scales from the + profile fit to the data. + Parameters ---------- model : DataModel @@ -287,19 +314,38 @@ def fit_profile(self, model): adjacent_inds = set(range(ind - self.pars['n_adjacent_cols'], ind + self.pars['n_adjacent_cols'] + 1)) adjacent_inds.discard(ind) valid_adjacent_inds = list(adjacent_inds.intersection(valid_profiles)) + # Cut out valid neighboring profiles - profile_data = model.data[self.custom_slice(dispaxis, valid_adjacent_inds)] + adjacent_condition = self.custom_slice(dispaxis, valid_adjacent_inds) + profile_data = model.data[adjacent_condition] + # Mask out bad pixels - profile_data = np.where( - model.dq[self.custom_slice(dispaxis, valid_adjacent_inds)] & self.DO_NOT_USE, - np.nan, - profile_data - ) - # Add additional cut to pull only from region with valid data for convenience (may not be necessary) - profile_data = profile_data[self.custom_slice(3 - dispaxis, list(range(profile_cut[0], profile_cut[1])))] + invalid_condition = (model.dq[adjacent_condition] & self.DO_NOT_USE).astype(bool) + profile_data[invalid_condition] = np.nan + + # Add additional cut to pull only from region with valid data + # for convenience (may not be necessary) + region_condition = self.custom_slice(3 - dispaxis, range(*profile_cut)) + profile_data = profile_data[region_condition] # Normalize profile data - normalized = profile_data / np.nanmax(np.abs(profile_data), axis=(dispaxis - 1), keepdims=True) + # TODO: check on signs here - absolute max sometimes picks up + # large negative outliers + profile_norm_scale = np.nanmax(np.abs(profile_data), axis=(dispaxis - 1), keepdims=True) + normalized = profile_data / profile_norm_scale + + # Get corresponding error and variance data and scale and mask to match + # Handle the variance arrays as errors, so the scales match. + err_names = ['err', 'var_poisson', 'var_rnoise', 'var_flat'] + norm_errors = {} + for err_name in err_names: + if err_name.startswith('var'): + err = np.sqrt(getattr(model, err_name)) + else: + err = getattr(model, err_name) + norm_err = err[adjacent_condition] + norm_err[invalid_condition] = np.nan + norm_errors[err_name] = norm_err[region_condition] / profile_norm_scale # Pull median for each pixel across profile. # Profile entry full of NaN values would produce a numpy @@ -307,12 +353,16 @@ def fit_profile(self, model): # so we suppress that above. median_profile = np.nanmedian(normalized, axis=(2 - dispaxis)) - # check_output[:, i] = median_profile + # Do the same for the errors + for err_name in norm_errors: + norm_errors[err_name] = np.nanmedian( + norm_errors[err_name], axis=(2 - dispaxis)) # Clean current profile of values flagged as bad - current_profile = model.data[self.custom_slice(dispaxis, ind)] + current_condition = self.custom_slice(dispaxis, ind) + current_profile = model.data[current_condition] cleaned_current = np.where( - model.dq[self.custom_slice(dispaxis, ind)] & self.DO_NOT_USE, + model.dq[current_condition] & self.DO_NOT_USE, np.nan, current_profile )[range(*profile_cut)] @@ -326,37 +376,76 @@ def fit_profile(self, model): norm_current = min_current / np.max(min_current) # Scale median profile to current profile with bad pixel - minimize mse? + # TODO: check on signs here - absolute max sometimes picks up + # large negative outliers norm_scale = minimize(self.profile_mse, x0=np.abs(np.nanmax(norm_current)), - args=(min_median, norm_current)).x + args=(min_median, norm_current)).x scale = np.max(min_current) + # Replace pixels that are do-not-use but not non-science + current_dq = model.dq[current_condition][range(*profile_cut)] + replace_condition = (current_dq & self.DO_NOT_USE + ^ current_dq & self.NON_SCIENCE) == 1 replaced_current = np.where( - (model.dq[self.custom_slice(dispaxis, ind)][range(*profile_cut)] & self.DO_NOT_USE ^ - model.dq[self.custom_slice(dispaxis, ind)][range(*profile_cut)] & self.NON_SCIENCE) == 1, + replace_condition, median_profile * norm_scale * scale, cleaned_current ) # Change the dq bits where old flag was DO_NOT_USE and new value is not nan - current_dq = model.dq[self.custom_slice(dispaxis, ind)][range(*profile_cut)] replaced_dq = np.where( - (current_dq & self.DO_NOT_USE ^ current_dq & self.NON_SCIENCE == 1) & - ~(np.isnan(replaced_current)), + replace_condition & ~(np.isnan(replaced_current)), current_dq ^ self.DO_NOT_USE ^ self.FLUX_ESTIMATED, current_dq ) - model_replaced.data[self.custom_slice(dispaxis, ind)][range(*profile_cut)] = replaced_current - model_replaced.dq[self.custom_slice(dispaxis, ind)][range(*profile_cut)] = replaced_dq + # Update data and DQ in the output model + model_replaced.data[current_condition][range(*profile_cut)] = replaced_current + model_replaced.dq[current_condition][range(*profile_cut)] = replaced_dq + + # Also update the errors and variances + current_err = model.err[current_condition][range(*profile_cut)] + replaced_err = np.where( + replace_condition, + norm_errors['err'] * norm_scale * scale, + current_err + ) + model_replaced.err[current_condition][range(*profile_cut)] = replaced_err + + current_var = model.var_poisson[current_condition][range(*profile_cut)] + replaced_var = np.where( + replace_condition, + (norm_errors['var_poisson'] * norm_scale * scale)**2, + current_var + ) + model_replaced.var_poisson[current_condition][range(*profile_cut)] = replaced_var + + current_var = model.var_rnoise[current_condition][range(*profile_cut)] + replaced_var = np.where( + replace_condition, + (norm_errors['var_rnoise'] * norm_scale * scale)**2, + current_var + ) + model_replaced.var_rnoise[current_condition][range(*profile_cut)] = replaced_var + + current_var = model.var_flat[current_condition][range(*profile_cut)] + replaced_var = np.where( + replace_condition, + (norm_errors['var_flat'] * norm_scale * scale)**2, + current_var + ) + model_replaced.var_flat[current_condition][range(*profile_cut)] = replaced_var return model_replaced def mingrad(self, model): """ - Minimum gradient replacement method; test the gradient along the spatial - and spectral axes using immediately adjacent pixels. Pick whichever dimension - has the minimum absolute gradient and replace the missing pixel with the average + Replace pixels with the minimum gradient replacement method. + + Test the gradient along the spatial and spectral axes using + immediately adjacent pixels. Pick whichever dimension has the minimum + absolute gradient and replace the missing pixel with the average of the two adjacent pixels along that dimension. This aims to make the process extremely local; near point sources it should do @@ -393,14 +482,22 @@ def mingrad(self, model): # Input data, err, and dq values indata = model.data - inerr = model.err indq = model.dq + inerr = model.err - # Output data, err, and dq values + # Propagate variance components as errors to get the scales right + in_var_p = np.sqrt(model.var_poisson) + in_var_r = np.sqrt(model.var_rnoise) + in_var_f = np.sqrt(model.var_flat) + + # Output data, err, var, and dq values model_replaced = model.copy() newdata = model_replaced.data - newerr = model_replaced.err newdq = model_replaced.dq + newerr = model_replaced.err + new_var_p = model_replaced.var_poisson + new_var_r = model_replaced.var_rnoise + new_var_f = model_replaced.var_flat # Make an array of x/y values on the detector (ysize, xsize) = indata.shape @@ -418,28 +515,50 @@ def mingrad(self, model): for ii in range(0, len(xindx)): left_data, right_data = indata[yindx[ii], xindx[ii] - 1], indata[yindx[ii], xindx[ii] + 1] top_data, bottom_data = indata[yindx[ii] - 1, xindx[ii]], indata[yindx[ii] + 1, xindx[ii]] + left_err, right_err = inerr[yindx[ii], xindx[ii] - 1], inerr[yindx[ii], xindx[ii] + 1] top_err, bottom_err = inerr[yindx[ii] - 1, xindx[ii]], inerr[yindx[ii] + 1, xindx[ii]] + left_var_p, right_var_p = in_var_p[yindx[ii], xindx[ii] - 1], in_var_p[yindx[ii], xindx[ii] + 1] + top_var_p, bottom_var_p = in_var_p[yindx[ii] - 1, xindx[ii]], in_var_p[yindx[ii] + 1, xindx[ii]] + + left_var_r, right_var_r = in_var_r[yindx[ii], xindx[ii] - 1], in_var_r[yindx[ii], xindx[ii] + 1] + top_var_r, bottom_var_r = in_var_r[yindx[ii] - 1, xindx[ii]], in_var_r[yindx[ii] + 1, xindx[ii]] + + left_var_f, right_var_f = in_var_f[yindx[ii], xindx[ii] - 1], in_var_f[yindx[ii], xindx[ii] + 1] + top_var_f, bottom_var_f = in_var_f[yindx[ii] - 1, xindx[ii]], in_var_f[yindx[ii] + 1, xindx[ii]] + # Compute absolute difference (slope) and average value in each direction (may be NaN) diffs = np.array([np.abs(left_data - right_data), np.abs(top_data - bottom_data)]) interp_data = np.array([(left_data + right_data) / 2., (top_data + bottom_data) / 2.]) interp_err = np.array([(left_err + right_err) / 2., (top_err + bottom_err) / 2.]) + interp_var_p = np.array([(left_var_p + right_var_p) / 2., (top_var_p + bottom_var_p) / 2.]) + interp_var_r = np.array([(left_var_r + right_var_r) / 2., (top_var_r + bottom_var_r) / 2.]) + interp_var_f = np.array([(left_var_f + right_var_f) / 2., (top_var_f + bottom_var_f) / 2.]) # Replace with the value from the lowest absolute slope estimator that was not NaN try: indmin = np.nanargmin(diffs) newdata[yindx[ii], xindx[ii]] = interp_data[indmin] newerr[yindx[ii], xindx[ii]] = interp_err[indmin] - # If original pixel was in the science array set DQ accordingly - if ((indq[yindx[ii], xindx[ii]] & self.NON_SCIENCE) == 0): - newdq[yindx[ii], xindx[ii]] = self.FLUX_ESTIMATED - # If original pixel was in non-science region set DQ accordingly - else: - newdq[yindx[ii], xindx[ii]] = self.FLUX_ESTIMATED + self.NON_SCIENCE + self.DO_NOT_USE + + # Square the interpolated errors back into variance + new_var_p[yindx[ii], xindx[ii]] = interp_var_p[indmin] ** 2 + new_var_r[yindx[ii], xindx[ii]] = interp_var_r[indmin] ** 2 + new_var_f[yindx[ii], xindx[ii]] = interp_var_f[indmin] ** 2 + + # If original pixel was in the science array, remove + # the DO_NOT_USE flag + if ((indq[yindx[ii], xindx[ii]] & self.DO_NOT_USE) + and not (indq[yindx[ii], xindx[ii]] & self.NON_SCIENCE)): + newdq[yindx[ii], xindx[ii]] -= self.DO_NOT_USE + + # Either way, add the FLUX_ESTIMATED flag + newdq[yindx[ii], xindx[ii]] |= self.FLUX_ESTIMATED nreplaced += 1 - except Exception: + + except (IndexError, ValueError): pass model_replaced.data = newdata @@ -459,7 +578,7 @@ def custom_slice(self, dispaxis, index): Using module-defined HORIZONTAL=1, VERTICAL=2 - index : int or slice + index : int or list Index or indices of cross-dispersion vectors to slice diff --git a/jwst/pixel_replace/tests/__init__.py b/jwst/pixel_replace/tests/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/jwst/pixel_replace/tests/test_pixel_replace.py b/jwst/pixel_replace/tests/test_pixel_replace.py new file mode 100644 index 0000000000..3ab91c7eec --- /dev/null +++ b/jwst/pixel_replace/tests/test_pixel_replace.py @@ -0,0 +1,240 @@ +import numpy as np +import pytest + +from stdatamodels.jwst import datamodels +from stdatamodels.jwst.datamodels.dqflags import pixel as flags + +from jwst.assign_wcs import AssignWcsStep +from jwst.assign_wcs.tests.test_nirspec import create_nirspec_ifu_file +from jwst.pixel_replace.pixel_replace_step import PixelReplaceStep + + +def cal_data(shape, bad_idx, dispaxis=1, model='slit'): + if model == 'image': + model = datamodels.ImageModel(shape) + elif model == 'ifu': + model = datamodels.IFUImageModel(shape) + else: + model = datamodels.SlitModel(shape) + model.meta.wcsinfo.dispersion_direction = dispaxis + + # Set the data and error arrays to all 1s except one bad pixel + # to correct at the middle of the array + model.data[:] = 1.0 + model.err[:] = 1.0 + model.var_poisson[:] = 1.0 + model.var_rnoise[:] = 1.0 + model.var_flat[:] = 1.0 + + bad_flag = flags['DO_NOT_USE'] + flags['OTHER_BAD_PIXEL'] + model.data[bad_idx] = np.nan + model.err[bad_idx] = np.nan + model.var_poisson[bad_idx] = np.nan + model.var_rnoise[bad_idx] = np.nan + model.var_flat[bad_idx] = np.nan + model.dq[bad_idx] = bad_flag + + # Also add a non-science region in one row and one column + non_science = flags['DO_NOT_USE']+ flags['NON_SCIENCE'] + model.data[..., 1] = np.nan + model.err[..., 1] = np.nan + model.var_poisson[..., 1] = np.nan + model.var_rnoise[..., 1] = np.nan + model.var_flat[..., 1] = np.nan + model.dq[..., 1] = non_science + + model.data[..., 1, :] = np.nan + model.err[..., 1, :] = np.nan + model.var_poisson[..., 1, :] = np.nan + model.var_rnoise[..., 1, :] = np.nan + model.var_flat[..., 1, :] = np.nan + model.dq[..., 1, :] = non_science + + return model + + +def nirspec_tso(): + bad_idx = (1, 10, 10) + model = cal_data(shape=(3, 20, 20), bad_idx=bad_idx, dispaxis=1) + model.meta.instrument.name = 'NIRSPEC' + model.meta.exposure.type = 'NRS_BRIGHTOBJ' + return model, bad_idx + + +def nirspec_fs_slitmodel(): + bad_idx = (10, 10) + model = cal_data(shape=(20, 20), bad_idx=bad_idx, dispaxis=1) + model.meta.instrument.name = 'NIRSPEC' + model.meta.exposure.type = 'NRS_FIXEDSLIT' + return model, bad_idx + + +def nirspec_msa_multislit(): + bad_idx = (10, 10) + slit_model = cal_data(shape=(20, 20), bad_idx=bad_idx, dispaxis=1) + model = datamodels.MultiSlitModel() + model.slits.append(slit_model) + model.meta.instrument.name = 'NIRSPEC' + model.meta.exposure.type = 'NRS_MSASPEC' + return model, bad_idx + + +def nirspec_ifu(): + shape = (2048, 2048) + bad_idx = (1414, 690) + + # IFU mode requires WCS information, so make a more realistic model + 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) + + model = datamodels.IFUImageModel(hdul) + model = AssignWcsStep.call(model) + + test_data = cal_data(shape=shape, bad_idx=bad_idx, dispaxis=1, model='ifu') + model.data = test_data.data + model.dq = test_data.dq + model.err = test_data.err + model.var_poisson = test_data.var_poisson + model.var_rnoise = test_data.var_rnoise + model.var_flat = test_data.var_flat + test_data.close() + + return model, bad_idx + + +def miri_lrs(): + bad_idx = (10, 10) + model = cal_data(shape=(20, 20), bad_idx=bad_idx, dispaxis=2, model='image') + model.meta.instrument.name = 'MIRI' + model.meta.exposure.type = 'MIR_LRS-FIXEDSLIT' + return model, bad_idx + + +def miri_mrs(): + shape = (20, 20) + + def mock_transform(*args): + return None, np.full(shape, 1), None + + bad_idx = (10, 10) + model = cal_data(shape=shape, bad_idx=bad_idx, dispaxis=2, model='ifu') + model.meta.instrument.name = 'MIRI' + model.meta.exposure.type = 'MIR_MRS' + model.meta.wcs = {'transform': mock_transform} + return model, bad_idx + + +@pytest.mark.parametrize('input_model_function', + [nirspec_tso, nirspec_fs_slitmodel, miri_lrs, miri_mrs]) +@pytest.mark.parametrize('algorithm', ['fit_profile', 'mingrad']) +def test_pixel_replace_no_container(input_model_function, algorithm): + """ + Test pixel replace for modes with no container. + + This includes ImageModel, SlitModel, and IFUImageModel. + """ + input_model, bad_idx = input_model_function() + + # for this simple case, the results from either algorithm should + # be the same + result = PixelReplaceStep.call(input_model, skip=False, algorithm=algorithm) + + for ext in ['data', 'err', 'var_poisson', 'var_rnoise', 'var_flat']: + # non-science edges are uncorrected + assert np.all(np.isnan(getattr(result, ext)[..., :, 1])) + assert np.all(np.isnan(getattr(result, ext)[..., 1, :])) + + # bad pixel is replaced: input had one nan value, output does not + assert np.isnan(getattr(input_model, ext)[bad_idx]) + assert getattr(result, ext)[bad_idx] == 1.0 + + # The DQ plane for the bad pixel is updated to remove do-not-use + # and add flux-estimated. The non-science edges are unchanged. + assert result.dq[bad_idx] == (input_model.dq[bad_idx] + - flags['DO_NOT_USE'] + + flags['FLUX_ESTIMATED']) + assert np.all(result.dq[..., :, 1] + == flags['DO_NOT_USE'] + flags['NON_SCIENCE']) + assert np.all(result.dq[..., 1, :] + == flags['DO_NOT_USE'] + flags['NON_SCIENCE']) + + result.close() + input_model.close() + + +@pytest.mark.parametrize('input_model_function', + [nirspec_msa_multislit]) +@pytest.mark.parametrize('algorithm', ['fit_profile', 'mingrad']) +def test_pixel_replace_multislit(input_model_function, algorithm): + """Test pixel replace for multislit modes.""" + input_model, bad_idx = input_model_function() + + # for this simple case, the results from either algorithm should + # be the same + result = PixelReplaceStep.call(input_model, skip=False, algorithm=algorithm) + + for ext in ['data', 'err', 'var_poisson', 'var_rnoise', 'var_flat']: + # non-science edges are uncorrected + assert np.all(np.isnan(getattr(result.slits[0], ext)[..., :, 1])) + assert np.all(np.isnan(getattr(result.slits[0], ext)[..., 1, :])) + + # bad pixel is replaced: input had one nan value, output does not + assert np.isnan(getattr(input_model.slits[0], ext)[bad_idx]) + assert getattr(result.slits[0], ext)[bad_idx] == 1.0 + + # The DQ plane for the bad pixel is updated to remove do-not-use + # and add flux-estimated. The non-science edges are unchanged. + assert result.slits[0].dq[bad_idx] == ( + input_model.slits[0].dq[bad_idx] + - flags['DO_NOT_USE'] + flags['FLUX_ESTIMATED']) + assert np.all(result.slits[0].dq[..., :, 1] + == flags['DO_NOT_USE'] + flags['NON_SCIENCE']) + assert np.all(result.slits[0].dq[..., 1, :] + == flags['DO_NOT_USE'] + flags['NON_SCIENCE']) + + result.close() + input_model.close() + + +@pytest.mark.slow +@pytest.mark.parametrize('input_model_function', + [nirspec_ifu]) +@pytest.mark.parametrize('algorithm', ['fit_profile', 'mingrad']) +def test_pixel_replace_nirspec_ifu(input_model_function, algorithm): + """ + Test pixel replacement for NIRSpec IFU. + + Larger data and more WCS operations required for testing make + this test take more than a minute, so marking this test 'slow'. + + The test is otherwise the same as for other modes. + """ + input_model, bad_idx = input_model_function() + + # for this simple case, the results from either algorithm should + # be the same + result = PixelReplaceStep.call(input_model, skip=False, algorithm=algorithm) + + for ext in ['data', 'err', 'var_poisson', 'var_rnoise', 'var_flat']: + # non-science edges are uncorrected + assert np.all(np.isnan(getattr(result, ext)[..., :, 1])) + assert np.all(np.isnan(getattr(result, ext)[..., 1, :])) + + # bad pixel is replaced: input had one nan value, output does not + assert np.isnan(getattr(input_model, ext)[bad_idx]) + assert getattr(result, ext)[bad_idx] == 1.0 + + # The DQ plane for the bad pixel is updated to remove do-not-use + # and add flux-estimated. The non-science edges are unchanged. + assert result.dq[bad_idx] == (input_model.dq[bad_idx] + - flags['DO_NOT_USE'] + + flags['FLUX_ESTIMATED']) + assert np.all(result.dq[..., :, 1] + == flags['DO_NOT_USE'] + flags['NON_SCIENCE']) + assert np.all(result.dq[..., 1, :] + == flags['DO_NOT_USE'] + flags['NON_SCIENCE']) + + result.close() + input_model.close() diff --git a/jwst/regtest/test_associations_sdp_pools.py b/jwst/regtest/test_associations_sdp_pools.py index 4f02c0692e..4444e40a3b 100644 --- a/jwst/regtest/test_associations_sdp_pools.py +++ b/jwst/regtest/test_associations_sdp_pools.py @@ -77,6 +77,11 @@ 'xfail': None, 'slow': True, }, + 'jw00818_20230407t030411_pool': { + 'args': [], + 'xfail': None, + 'slow': False, + }, 'jw00839_20221220t025418_pool': { 'args': ['-i', 'o002', 'c1000'], 'xfail': None, diff --git a/jwst/resample/resample_spec_step.py b/jwst/resample/resample_spec_step.py index fa705f5372..e215a7bb61 100755 --- a/jwst/resample/resample_spec_step.py +++ b/jwst/resample/resample_spec_step.py @@ -62,35 +62,12 @@ def process(self, input): output = input_new.meta.filename self.blendheaders = False - # Get the drizpars reference file - for reftype in self.reference_file_types: - ref_filename = self.get_reference_file(input_models[0], reftype) - - if ref_filename != 'N/A': - self.log.info('Drizpars reference file: {}'.format(ref_filename)) - kwargs = self.get_drizpars(ref_filename, input_models) - else: - # Deal with NIRSpec, which currently has no default drizpars reffile - self.log.info("No DRIZPARS reffile") - kwargs = self._set_spec_defaults() - kwargs['blendheaders'] = self.blendheaders - - kwargs['output_shape'] = self._check_list_pars( - self.output_shape, - 'output_shape', - min_vals=[1, 1] - ) - kwargs['output_wcs'] = self._load_custom_wcs( - self.output_wcs, - kwargs['output_shape'] - ) - - kwargs['allowed_memory'] = self.allowed_memory + # Setup drizzle-related parameters + kwargs = self.get_drizpars() kwargs['output'] = output - - # Call resampling self.drizpars = kwargs + # Call resampling if isinstance(input_models[0], MultiSlitModel): result = self._process_multislit(input_models) diff --git a/jwst/resample/resample_step.py b/jwst/resample/resample_step.py index f19e3cfcb9..1f99d12a41 100755 --- a/jwst/resample/resample_step.py +++ b/jwst/resample/resample_step.py @@ -2,10 +2,7 @@ import re from copy import deepcopy -import numpy as np import asdf -from astropy.extern.configobj.validate import Validator -from astropy.extern.configobj.configobj import ConfigObj from stdatamodels.jwst import datamodels @@ -43,24 +40,24 @@ class ResampleStep(Step): class_alias = "resample" spec = """ - pixfrac = float(default=1.0) # change back to None when drizpar reference files are updated - kernel = string(default='square') # change back to None when drizpar reference files are updated - fillval = string(default='NAN' ) - weight_type = option('ivm', 'exptime', None, default='ivm') # change back to None when drizpar ref update + pixfrac = float(min=0.0, max=1.0, default=1.0) # Pixel shrinkage factor + kernel = option('square','gaussian','point','turbo','lanczos2','lanczos3',default='square') # Flux distribution kernel + fillval = string(default='NAN') # Output value for pixels with no weight or flux + weight_type = option('ivm', 'exptime', None, default='ivm') # Input image weighting type output_shape = int_list(min=2, max=2, default=None) # [x, y] order crpix = float_list(min=2, max=2, default=None) crval = float_list(min=2, max=2, default=None) - rotation = float(default=None) - pixel_scale_ratio = float(default=1.0) # Ratio of input to output pixel scale - pixel_scale = float(default=None) # Absolute pixel scale in arcsec - output_wcs = string(default='') # Custom output WCS. - single = boolean(default=False) - blendheaders = boolean(default=True) - allowed_memory = float(default=None) # Fraction of memory to use for the combined image. - in_memory = boolean(default=True) + rotation = float(default=None) # Output image Y-axis PA relative to North + pixel_scale_ratio = float(default=1.0) # Ratio of input to output pixel scale + pixel_scale = float(default=None) # Absolute pixel scale in arcsec + output_wcs = string(default='') # Custom output WCS + single = boolean(default=False) # Resample each input to its own output grid + blendheaders = boolean(default=True) # Blend metadata from inputs into output + allowed_memory = float(default=None) # Fraction of memory to use for the combined image + in_memory = boolean(default=True) # Keep images in memory """ - reference_file_types = ['drizpars'] + reference_file_types = [] def process(self, input): @@ -87,39 +84,8 @@ def process(self, input): # resample can only handle 2D images, not 3D cubes, etc raise RuntimeError("Input {} is not a 2D image.".format(input_models[0])) - # Get drizzle parameters reference file, if there is one - self.wht_type = self.weight_type - if 'drizpars' in self.reference_file_types: - ref_filename = self.get_reference_file(input_models[0], 'drizpars') - else: # no drizpars reference file found - ref_filename = 'N/A' - - if ref_filename == 'N/A': - self.log.info('No drizpars reference file found.') - kwargs = self._set_spec_defaults() - else: - self.log.info('Using drizpars reference file: {}'.format(ref_filename)) - kwargs = self.get_drizpars(ref_filename, input_models) - - kwargs['allowed_memory'] = self.allowed_memory - - # Custom output WCS parameters. - # Modify get_drizpars if any of these get into reference files: - kwargs['output_shape'] = self._check_list_pars( - self.output_shape, - 'output_shape', - min_vals=[1, 1] - ) - kwargs['output_wcs'] = self._load_custom_wcs( - self.output_wcs, - kwargs['output_shape'] - ) - kwargs['crpix'] = self._check_list_pars(self.crpix, 'crpix') - kwargs['crval'] = self._check_list_pars(self.crval, 'crval') - kwargs['rotation'] = self.rotation - kwargs['pscale'] = self.pixel_scale - kwargs['pscale_ratio'] = self.pixel_scale_ratio - kwargs['in_memory'] = self.in_memory + # Setup drizzle-related parameters + kwargs = self.get_drizpars() # Call the resampling routine resamp = resample.ResampleData(input_models, output=output, **kwargs) @@ -195,128 +161,42 @@ def _load_custom_wcs(asdf_wcs_file, output_shape): return wcs - def get_drizpars(self, ref_filename, input_models): + def get_drizpars(self): """ - Extract drizzle parameters from reference file. - - This method extracts parameters from the drizpars reference file and - uses those to set defaults on the following ResampleStep configuration - parameters: - - pixfrac = float(default=None) - kernel = string(default=None) - fillval = string(default=None) - wht_type = option('ivm', 'exptime', None, default=None) - - Once the defaults are set from the reference file, if the user has - used a resample.cfg file or run ResampleStep using command line args, - then these will overwrite the defaults pulled from the reference file. + Load all drizzle-related parameter values into kwargs list. """ - with datamodels.DrizParsModel(ref_filename) as drpt: - drizpars_table = drpt.data - - num_groups = len(input_models.group_names) - filtname = input_models[0].meta.instrument.filter - row = None - filter_match = False - # look for row that applies to this set of input data models - for n, filt, num in zip( - range(0, len(drizpars_table)), - drizpars_table['filter'], - drizpars_table['numimages'] - ): - # only remember this row if no exact match has already been made for - # the filter. This allows the wild-card row to be anywhere in the - # table; since it may be placed at beginning or end of table. - - if str(filt) == "ANY" and not filter_match and num_groups >= num: - row = n - # always go for an exact match if present, though... - if filtname == filt and num_groups >= num: - row = n - filter_match = True - - # With presence of wild-card rows, code should never trigger this logic - if row is None: - self.log.error("No row found in %s matching input data.", ref_filename) - raise ValueError - - # Define the keys to pull from drizpars reffile table. - # All values should be None unless the user set them on the command - # line or in the call to the step - - drizpars = dict( + # Define the keys pulled from step parameters + kwargs = dict( pixfrac=self.pixfrac, kernel=self.kernel, fillval=self.fillval, - wht_type=self.weight_type - # pscale_ratio=self.pixel_scale_ratio, # I think this can be removed JEM (??) - ) - - # For parameters that are set in drizpars table but not set by the - # user, use these. Otherwise, use values set by user. - reffile_drizpars = {k: v for k, v in drizpars.items() if v is None} - user_drizpars = {k: v for k, v in drizpars.items() if v is not None} - - # read in values from that row for each parameter - for k in reffile_drizpars: - if k in drizpars_table.names: - reffile_drizpars[k] = drizpars_table[k][row] - - # Convert the strings in the FITS binary table from np.bytes_ to str - for k, v in reffile_drizpars.items(): - if isinstance(v, np.bytes_): - reffile_drizpars[k] = v.decode('UTF-8') - - all_drizpars = {**reffile_drizpars, **user_drizpars} - - kwargs = dict( + wht_type=self.weight_type, good_bits=GOOD_BITS, single=self.single, - blendheaders=self.blendheaders + blendheaders=self.blendheaders, + allowed_memory = self.allowed_memory, + in_memory = self.in_memory ) - kwargs.update(all_drizpars) - - for k, v in kwargs.items(): - self.log.debug(' {}={}'.format(k, v)) - - return kwargs - - def _set_spec_defaults(self): - """NIRSpec currently has no default drizpars reference file, so default - drizzle parameters are not set properly. This method sets them. - - Remove this class method when a drizpars reffile is delivered. - """ - configspec = self.load_spec_file() - config = ConfigObj(configspec=configspec) - if config.validate(Validator()): - kwargs = config.dict() - - if self.pixfrac is None: - self.pixfrac = 1.0 - if self.kernel is None: - self.kernel = 'square' - if self.fillval is None: - self.fillval = 'NAN' - # Force definition of good bits - kwargs['good_bits'] = GOOD_BITS - kwargs['pixfrac'] = self.pixfrac - kwargs['kernel'] = str(self.kernel) - kwargs['fillval'] = str(self.fillval) - # self.weight_type has a default value of None - # The other instruments read this parameter from a reference file - if self.wht_type is None: - self.wht_type = 'ivm' - - kwargs['wht_type'] = str(self.wht_type) + # Custom output WCS parameters. + kwargs['output_shape'] = self._check_list_pars( + self.output_shape, + 'output_shape', + min_vals=[1, 1] + ) + kwargs['output_wcs'] = self._load_custom_wcs( + self.output_wcs, + kwargs['output_shape'] + ) + kwargs['crpix'] = self._check_list_pars(self.crpix, 'crpix') + kwargs['crval'] = self._check_list_pars(self.crval, 'crval') + kwargs['rotation'] = self.rotation + kwargs['pscale'] = self.pixel_scale kwargs['pscale_ratio'] = self.pixel_scale_ratio - kwargs.pop('pixel_scale_ratio') + # Report values to processing log for k, v in kwargs.items(): - if k in ['pixfrac', 'kernel', 'fillval', 'wht_type', 'pscale_ratio']: - log.info(' using: %s=%s', k, repr(v)) + self.log.debug(' {}={}'.format(k, v)) return kwargs diff --git a/jwst/wavecorr/tests/test_wavecorr.py b/jwst/wavecorr/tests/test_wavecorr.py index e7681fb6d2..85e74a0d0a 100644 --- a/jwst/wavecorr/tests/test_wavecorr.py +++ b/jwst/wavecorr/tests/test_wavecorr.py @@ -7,6 +7,7 @@ from stdatamodels.jwst import datamodels from stdatamodels.jwst.transforms import models +from stpipe.crds_client import reference_uri_to_cache_path import jwst from jwst.assign_wcs import AssignWcsStep @@ -59,7 +60,7 @@ def test_wavecorr(): # test the round-tripping on the wavelength correction transform ref_name = im_wave.meta.ref_file.wavecorr.name freference = datamodels.WaveCorrModel( - WavecorrStep.reference_uri_to_cache_path(ref_name, im.crds_observatory)) + reference_uri_to_cache_path(ref_name, im.crds_observatory)) lam_uncorr = lam_before * 1e-6 wave2wavecorr = wavecorr.calculate_wavelength_correction_transform( @@ -160,7 +161,7 @@ def test_skipped(): ref_name = outw.meta.ref_file.wavecorr.name reffile = datamodels.WaveCorrModel( - WavecorrStep.reference_uri_to_cache_path(ref_name, im.crds_observatory)) + reference_uri_to_cache_path(ref_name, im.crds_observatory)) source_xpos = 0.1 aperture_name = 'S200A1' @@ -263,10 +264,10 @@ def test_wavecorr_fs(): # test the round-tripping on the wavelength correction transform ref_name = result.meta.ref_file.wavecorr.name - freference = datamodels.WaveCorrModel(WavecorrStep.reference_uri_to_cache_path(ref_name, im.crds_observatory)) + freference = datamodels.WaveCorrModel(reference_uri_to_cache_path(ref_name, im.crds_observatory)) lam_uncorr = lam_before * 1e-6 - wave2wavecorr = wavecorr.calculate_wavelength_correction_transform(lam_uncorr, dispersion, + wave2wavecorr = wavecorr.calculate_wavelength_correction_transform(lam_uncorr, dispersion, freference, slit.source_xpos, 'S200A1') lam_corr = wave2wavecorr(lam_uncorr) assert_allclose(lam_uncorr, wave2wavecorr.inverse(lam_corr))