diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index f9f78fea0d..cb1591fdeb 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -21,7 +21,7 @@ jobs: build: needs: [ check ] if: always() && (needs.check.result == 'success' || needs.check.result == 'skipped') - uses: OpenAstronomy/github-actions-workflows/.github/workflows/publish.yml@924441154cf3053034c6513d5e06c69d262fb9a6 # v1.13.0 + uses: OpenAstronomy/github-actions-workflows/.github/workflows/publish.yml@9f1f43251dde69da8613ea8e11144f05cdea41d5 # v1.15.0 with: upload_to_pypi: ${{ (github.event_name == 'release') && (github.event.action == 'released') }} targets: | diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index daf0c4221b..1d712bc1a6 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -40,14 +40,15 @@ jobs: python-version: '3.12' - uses: pre-commit/action@2c7b3805fd2a0fd8c1884dcaebf91fc102a13ecd # v3.0.1 check: - uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@924441154cf3053034c6513d5e06c69d262fb9a6 # v1.13.0 + uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@9f1f43251dde69da8613ea8e11144f05cdea41d5 # v1.15.0 with: default_python: "3.12" envs: | - linux: check-dependencies - linux: check-types test: - uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@924441154cf3053034c6513d5e06c69d262fb9a6 # v1.13.0 + uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@9f1f43251dde69da8613ea8e11144f05cdea41d5 # v1.15.0 + needs: [ crds_context ] with: setenv: | CRDS_PATH: /tmp/data/crds_cache diff --git a/.github/workflows/ci_cron.yml b/.github/workflows/ci_cron.yml index 535a93d878..3f19dbe425 100644 --- a/.github/workflows/ci_cron.yml +++ b/.github/workflows/ci_cron.yml @@ -24,7 +24,8 @@ concurrency: jobs: test: if: (github.repository == 'spacetelescope/jwst' && (github.event_name == 'schedule' || github.event_name == 'push' || github.event_name == 'workflow_dispatch' || contains(github.event.pull_request.labels.*.name, 'run scheduled tests'))) - uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@924441154cf3053034c6513d5e06c69d262fb9a6 # v1.13.0 + uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@9f1f43251dde69da8613ea8e11144f05cdea41d5 # v1.15.0 + needs: [ crds_context ] with: setenv: | CRDS_PATH: /tmp/crds_cache diff --git a/.github/workflows/tests_devdeps.yml b/.github/workflows/tests_devdeps.yml index 446b74f1b9..8db7c09c01 100644 --- a/.github/workflows/tests_devdeps.yml +++ b/.github/workflows/tests_devdeps.yml @@ -33,7 +33,8 @@ concurrency: jobs: test: if: (github.repository == 'spacetelescope/jwst' && (github.event_name == 'schedule' || github.event_name == 'push' || github.event_name == 'workflow_dispatch' || contains(github.event.pull_request.labels.*.name, 'run devdeps tests'))) - uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@924441154cf3053034c6513d5e06c69d262fb9a6 # v1.13.0 + uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@9f1f43251dde69da8613ea8e11144f05cdea41d5 # v1.15.0 + needs: [ crds_context ] with: setenv: | CRDS_PATH: /tmp/data/crds_cache diff --git a/changes/8846.ami.rst b/changes/8846.ami.rst new file mode 100644 index 0000000000..a965c823ff --- /dev/null +++ b/changes/8846.ami.rst @@ -0,0 +1 @@ +Change how AMI observables are averaged: average fringe quantities before calculating additional observables. Update their error calculation: use covariance of amplitudes/phases (and derived quantities) and standard error of the mean. Code now expects an ASDF filename string for user-supplied affine2d and bandpass input arguments. Example file creation in documentation. \ No newline at end of file diff --git a/changes/8847.master_background.rst b/changes/8847.master_background.rst new file mode 100644 index 0000000000..f1b516a1f2 --- /dev/null +++ b/changes/8847.master_background.rst @@ -0,0 +1 @@ +Include resample, pixel_replace and extract_1d into MOS background pipeline diff --git a/changes/8890.jump.rst b/changes/8890.jump.rst new file mode 100644 index 0000000000..d04413a547 --- /dev/null +++ b/changes/8890.jump.rst @@ -0,0 +1,2 @@ +Add maximum_shower_amplitude parameter to MIRI cosmic rays showers routine +to fix accidental flagging of bright science pixels. diff --git a/changes/8957.general.rst b/changes/8957.general.rst new file mode 100644 index 0000000000..91a475dc0f --- /dev/null +++ b/changes/8957.general.rst @@ -0,0 +1 @@ +remove ``okify_regtests`` script (move to ``ci_watson``) diff --git a/changes/8966.general.rst b/changes/8966.general.rst new file mode 100644 index 0000000000..5b21ca6897 --- /dev/null +++ b/changes/8966.general.rst @@ -0,0 +1 @@ +Increase asdf upper pin to 5 diff --git a/changes/8978.assign_wcs.rst b/changes/8978.assign_wcs.rst new file mode 100644 index 0000000000..2771c1d3c6 --- /dev/null +++ b/changes/8978.assign_wcs.rst @@ -0,0 +1 @@ +Fix negative SHUTTRID values under numpy 2 diff --git a/docs/jwst/ami_analyze/description.rst b/docs/jwst/ami_analyze/description.rst index c0e2d27016..26b2f597ff 100644 --- a/docs/jwst/ami_analyze/description.rst +++ b/docs/jwst/ami_analyze/description.rst @@ -36,7 +36,7 @@ other options: rotation search values. The default setting of '-3 3 1' results in search values of [-3, -2, -1, 0, 1, 2, 3]. -:--bandpass: Synphot spectrum or suitable array to override filter/source +:--bandpass: ASDF file containing suitable array to override filter/source (default=None) :--usebp: If True, exclude pixels marked DO_NOT_USE from fringe fitting @@ -47,10 +47,73 @@ other options: :--chooseholes: If not None, fit only certain fringes e.g. ['B4','B5','B6','C2'] (default=None) -:--affine2d: User-defined Affine2d object (default=None) +:--affine2d: ASDF file containing user-defined affine parameters (default=None) :--run_bpfix: Run Fourier bad pixel fix on cropped data (default=True) - + + +Creating ASDF files +^^^^^^^^^^^^^^^^^^^ +The optional arguments `bandpass` and `affine2d` must be written to `ASDF `_ +files to be used by the step. The step expects the contents to be stored with particular keys but the format is not currently +enforced by a schema; incorrect ASDF file contents will cause the step to revert back to the defaults for each argument. + +Examples of how to create ASDF files containing the properly formatted information for each of the arguments follows. + +.. code-block:: python + + # Create a F480M filter + Vega bandpass ASDF file + + import asdf + from jwst.ami import utils + from stdatamodels.jwst import datamodels + from synphot import SourceSpectrum + + # F480M throughput reference file from JWST CRDS + throughput_file = 'jwst_niriss_throughput_0012.fits' + nspecbin=19 + throughput_model = datamodels.open(throughput_file) + + filt_spec = utils.get_filt_spec(throughput_model) + src_spec = SourceSpectrum.from_vega() + bandpass = utils.combine_src_filt(filt_spec, + src_spec, + trim=0.01, + nlambda=nspecbin) + + # This bandpass has shape (19, 2); each row is [throughput, wavelength] + asdf_name = 'bandpass_f480m_vega.asdf' + tree = {"bandpass": bandpass} + with open(asdf_name, 'wb') as fh: + af = asdf.AsdfFile(tree) + af.write_to(fh) + af.close() + throughput_model.close() + + +.. code-block:: python + + # Create an affine transform ASDF file to use for the model + + import asdf + tree = { + 'mx': 1., # dimensionless x-magnification + 'my': 1., # dimensionless y-magnification + 'sx': 0., # dimensionless x shear + 'sy': 0., # dimensionless y shear + 'xo': 0., # x-offset in pupil space + 'yo': 0., # y-offset in pupil space + 'rotradccw': None + } + + affineasdf = 'affine.asdf' + + with open(affineasdf, 'wb') as fh: + af = asdf.AsdfFile(tree) + af.write_to(fh) + af.close() + + Inputs ------ @@ -70,7 +133,7 @@ The ``ami_analyze`` step itself does not accept an ASN as input. Outputs ------- -The ``ami_analyze`` step produces three output files. The first two (``_ami-oi.fits`` and ``_amimulti-oi.fits``) contain the interferometric observables, and the third (``_amilg.fits``) contains the data, LG model, and residuals. These are described in more detail below. +The ``ami_analyze`` step produces three output files. The first two (``_ami-oi.fits`` and ``_amimulti-oi.fits``) contain the interferometric observables, and the third (``_amilg.fits``) contains the data, LG model, and residuals. All products contain a ``PRIMARY`` HDU containing header keywords but no science data, as well as an ``ASDF`` extension. The files are described in more detail below. The output file name syntax is exposure-based, using the input file name as the root, with the addition of the association candidate ID and the "_ami-oi", "_amimulti-oi", or "amilg" product type suffix, e.g. @@ -81,20 +144,19 @@ Interferometric observables :Data model: `~jwst.datamodels.AmiOIModel` :File suffix: _ami-oi.fits, _amimulti-oi.fits -. The inteferometric observables are saved as OIFITS files, a registered FITS format +The inteferometric observables are saved as OIFITS files, a registered FITS format for optical interferometry, containing the following list of extensions: 1) ``OI_ARRAY``: AMI subaperture information 2) ``OI_TARGET``: target properties -3) ``OI_T3``: extracted closure amplitudes, phases +3) ``OI_T3``: extracted closure amplitudes, triple-product phases 4) ``OI_VIS``: extracted visibility (fringe) amplitudes, phases 5) ``OI_VIS2``: squared visibility (fringe) amplitudes 6) ``OI_WAVELENGTH``: filter information For more information on the format and contents of OIFITS files, see the `OIFITS2 standard `_. -The _ami-oi.fits file contains tables of median observables over all integrations of the input file. Errors -are computed as the sigma-clipped standard deviation over integrations. +The _ami-oi.fits file contains tables of observables averaged over all integrations of the input file. The error is taken to be the standard error of the mean, where the variance is the covariance between amplitudes and phases (e.g. fringe amplitudes and fringe phases, closure phases and triple-product amplitudes). The _amimulti-oi.fits file contains observables for each integration, and does not contain error estimates. The structure is the same as the _ami-oi.fits file, but the following data columns are 2D, with the second dimension being the number of integrations: "PISTONS", "PIST_ERR", "VISAMP", "VISAMPERR", "VISPHI", "VISPHIERR", "VIS2DATA", "VIS2ERR", "T3AMP", "T3AMPERR", "T3PHI", "T3PHIERR". diff --git a/docs/jwst/jump/arguments.rst b/docs/jwst/jump/arguments.rst index 9eada8a546..778d42363b 100644 --- a/docs/jwst/jump/arguments.rst +++ b/docs/jwst/jump/arguments.rst @@ -111,6 +111,8 @@ is defined as: * ``--min_diffs_single_pass``: The minimum number of differences to switch to flagging all outliers at once +* ``--max_shower_amplitude``: The maximum allowable amplitude for MIRI showers in DN/s + **Parameter that affects both Snowball and Shower flagging** * ``--max_extended_radius``: The maxiumum extension of the jump and saturation that will be flagged for showers or snowballs diff --git a/docs/jwst/jump/description.rst b/docs/jwst/jump/description.rst index c79f46c141..e7ac9aac9e 100644 --- a/docs/jwst/jump/description.rst +++ b/docs/jwst/jump/description.rst @@ -49,7 +49,8 @@ detectors are almost always circles with a central region that is saturated. The saturated core allows the search for smaller events without false positives. The mid-IR (MIRI) detectors do not, in general, have a saturated center and are only rarely circular. Thus, we fit the minimum enclosing ellipse and do not require that there are saturated pixels -within the ellipse. +within the ellipse. Likewise, MIRI showers are only flagged when detected features are consistent +with the maximum known amplitude (in DN/s) of shower artifacts. Multiprocessing --------------- diff --git a/docs/jwst/master_background/description.rst b/docs/jwst/master_background/description.rst index 96f90e82aa..d54b1eb0e6 100644 --- a/docs/jwst/master_background/description.rst +++ b/docs/jwst/master_background/description.rst @@ -376,6 +376,7 @@ as follows: extended sources (appropriate for background signal), and saving the extended source correction arrays for each slit in an internal copy of the data model #. If a user-supplied master background spectrum is **not** given, the + :ref:`pixel_replace `, :ref:`resample_spec ` and :ref:`extract_1d ` steps are applied to the calibrated background slits, resulting in extracted 1D background spectra diff --git a/docs/jwst/pipeline/calwebb_ami3.rst b/docs/jwst/pipeline/calwebb_ami3.rst index c3f494861f..d8aeb69cbc 100644 --- a/docs/jwst/pipeline/calwebb_ami3.rst +++ b/docs/jwst/pipeline/calwebb_ami3.rst @@ -7,8 +7,7 @@ calwebb_ami3: Stage 3 Aperture Masking Interferometry (AMI) Processing :Alias: calwebb_ami3 The stage 3 AMI pipeline is applied to associations of calibrated NIRISS AMI exposures. -It computes fringe parameters for individual exposures, averages the fringe results from -multiple exposures, and, optionally, corrects science target fringe parameters using the +It computes fringe parameters for individual exposures, and, optionally, corrects science target fringe parameters using the fringe results from reference PSF targets. The steps applied by the ``calwebb_ami3`` pipeline are shown below. @@ -83,15 +82,14 @@ Interferometric observables :Data model: `~jwst.datamodels.AmiOIModel` :File suffix: _ami-oi.fits - -For every input exposure, the fringe parameters and closure phases calculated +For every input exposure, the interferometric observables calculated by the :ref:`ami_analyze ` step are saved to an "_ami-oi.fits" product file, -which is a FITS table of median observables over all integrations of the input file. +which is a FITS table of averaged observables over all integrations of the input file. Product names use the input "_calints" exposure-based file name, with the association candidate ID included and the product type changed to "_ami-oi.fits", e.g. "jw93210001001_03101_00001_nis_a0003_ami-oi.fits." -Normalized Interferometric Observables +Normalized interferometric observables ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ :Data model: `~jwst.datamodels.AmiOIModel` :File suffix: _aminorm-oi.fits @@ -102,3 +100,8 @@ via the :ref:`ami_normalize ` step, and will be saved to an product file. This file has the same FITS table format as the "_ami-oi.fits" products. The file name root uses the source-based output product name given in the ASN file, e.g. "jw93210-a0003_t001_niriss_f480m-nrm_aminorm-oi.fits." + +.. note:: + + Users may wish to run the :ref:`ami_analyze step ` separately for access to flexible input arguments and to save additional diagnostic output products. See the step documentation for more details. + diff --git a/jwst/ami/ami_analyze_step.py b/jwst/ami/ami_analyze_step.py index ddf6d3fd5f..3330e6aa25 100755 --- a/jwst/ami/ami_analyze_step.py +++ b/jwst/ami/ami_analyze_step.py @@ -2,6 +2,10 @@ from ..stpipe import Step from . import ami_analyze +from . import utils + +import numpy as np +import asdf __all__ = ["AmiAnalyzeStep"] @@ -16,11 +20,11 @@ class AmiAnalyzeStep(Step): rotation = float(default=0.0) # Rotation initial guess [deg] psf_offset = string(default='0.0 0.0') # PSF offset values to use to create the model array rotation_search = string(default='-3 3 1') # Rotation search parameters: start, stop, step - bandpass = any(default=None) # Synphot spectrum or array to override filter/source + bandpass = string(default=None) # ASDF file containing array to override filter/source usebp = boolean(default=True) # If True, exclude pixels marked DO_NOT_USE from fringe fitting firstfew = integer(default=None) # If not None, process only the first few integrations chooseholes = string(default=None) # If not None, fit only certain fringes e.g. ['B4','B5','B6','C2'] - affine2d = any(default=None) # None or user-defined Affine2d object + affine2d = string(default=None) # ASDF file containing user-defined affine parameters run_bpfix = boolean(default=True) # Run Fourier bad pixel fix on cropped data """ @@ -32,6 +36,96 @@ def save_model(self, model, *args, **kwargs): kwargs["suffix"] = ["ami-oi", "amimulti-oi", "amilg"][kwargs.pop("idx")] return Step.save_model(self, model, *args, **kwargs) + def override_bandpass(self): + """ + Read bandpass from asdf file and use it to override the default. + + Expects an array of [effstims, wave_m] + (i.e. np.array((effstims,wave_m)).T) stored as 'bandpass' in asdf file, + where effstims are normalized countrates (unitless) and wave_m are the + wavelengths across the filter at which to compute the model (meters). + + Returns + ------- + bandpass: array + Array of [countrates, wavelengths] + """ + + try: + with asdf.open(self.bandpass, lazy_load=False) as af: + bandpass = np.array(af['bandpass']) + + # assume it is an array of the correct shape + wavemin = np.min(bandpass[:,1]) + wavemax = np.max(bandpass[:,1]) + self.log.info('User-defined bandpass provided:') + self.log.info('\tOVERWRITING ALL NIRISS-SPECIFIC FILTER/BANDPASS VARIABLES') + self.log.info(f'Using {bandpass.shape[0]} wavelengths for fit.') + self.log.info(f'Wavelength min: {wavemin:.3e} \t Wavelength max: {wavemax:.3e}') + + # update attribute and return + self.bandpass = bandpass + return bandpass + + except FileNotFoundError: + message = f'File {self.bandpass} could not be found at the specified location.' + raise Exception(message) + + except KeyError: + message1 = 'ASDF file does not contain the required "bandpass" key. ' + message2 = 'See step documentation for info on creating a custom bandpass ASDF file.' + raise Exception((message1 + message2)) + + except (IndexError, ValueError): + message1 = f'Could not use bandpass from {self.bandpass}. It may have the wrong shape. ' + message2 = 'See documentation for info on creating a custom bandpass ASDF file.' + raise Exception((message1 + message2)) + + def override_affine2d(self): + """ + Read user-input affine transform from ASDF file. + + Makes an Affine2d object (see utils.Affine2D class). + Input should contain mx,my,sx,sy,xo,yo,rotradccw. + """ + try: + with asdf.open(self.affine2d, lazy_load=False) as af: + affine2d = utils.Affine2d( + mx = af['mx'], + my = af['my'], + sx = af['sx'], + sy = af['sy'], + xo = af['xo'], + yo = af['yo'], + rotradccw = af['rotradccw'] + ) + self.log.info(f'Using affine transform from ASDF file {self.affine2d}') + # now self.affine2d updated from string to object + self.affine2d = affine2d + return affine2d + + except FileNotFoundError: + self.log.info(f'File {self.affine2d} could not be found at the specified location.') + self.log.info('\t **** DEFAULTING TO USE IDENTITY TRANSFORM ****') + affine2d = None + + except KeyError: + message1 = 'ASDF file does not contain all of the required keys: mx, my, sx, sy ,xo, yo, rotradccw. ' + message2 = 'See step documentation for info on creating a custom affine2d ASDF file.' + self.log.info((message1 + message2)) + self.log.info('\t **** DEFAULTING TO USE IDENTITY TRANSFORM ****') + affine2d = None + + except (IndexError, TypeError, ValueError): + message1 = f'Could not use affine2d from {self.affine2d}. ' + message2 = 'See documentation for info on creating a custom bandpass ASDF file.' + self.log.info((message1 + message2)) + self.log.info('\t **** DEFAULTING TO USE IDENTITY TRANSFORM ****') + affine2d = None + + self.affine2d = affine2d + return affine2d + def process(self, input): """ Performs analysis of an AMI mode exposure by applying the LG algorithm. @@ -83,6 +177,15 @@ def process(self, input): raise RuntimeError("No THROUGHPUT reference file found. " "ami_analyze cannot continue.") + # If there's a user-defined bandpass or affine, handle it + if bandpass is not None: + bandpass = self.override_bandpass() + + if affine2d is not None: + # if it is None, handled in apply_LG_plus + affine2d = self.override_affine2d() + + # Get the name of the NRM reference file to use nrm_reffile = self.get_reference_file(input_model, 'nrm') self.log.info(f'Using NRM reference file {nrm_reffile}') diff --git a/jwst/ami/instrument_data.py b/jwst/ami/instrument_data.py index 7b5e11614d..26d4f73d08 100755 --- a/jwst/ami/instrument_data.py +++ b/jwst/ami/instrument_data.py @@ -179,10 +179,11 @@ def read_data_model(self, input_model): pscale_deg = np.mean([pscaledegx, pscaledegy]) self.pscale_rad = np.deg2rad(pscale_deg) self.pscale_mas = pscale_deg * (60 * 60 * 1000) - self.pav3 = input_model.meta.pointing.pa_v3 + + self.roll_ref = input_model.meta.wcsinfo.roll_ref self.vparity = input_model.meta.wcsinfo.vparity self.v3iyang = input_model.meta.wcsinfo.v3yangle - self.parangh = input_model.meta.wcsinfo.roll_ref + self.crpix1 = input_model.meta.wcsinfo.crpix1 self.crpix2 = input_model.meta.wcsinfo.crpix2 self.pupil = input_model.meta.instrument.pupil @@ -223,7 +224,7 @@ def read_data_model(self, input_model): log.warning("All integrations will be analyzed") self.nwav = scidata.shape[0] [self.wls.append(self.wls[0]) for f in range(self.nwav - 1)] - # Rotate mask hole centers by pav3 + v3i_yang to be in equatorial coordinates + # Rotate mask hole centers by roll_ref + v3i_yang to be in equatorial coordinates ctrs_sky = self.mast2sky() oifctrs = np.zeros(self.mask.ctrs.shape) oifctrs[:, 0] = ctrs_sky[:, 1].copy() * -1 @@ -331,11 +332,19 @@ def reset_nwav(self, nwav): def mast2sky(self): """ - Rotate hole center coordinates: - Clockwise by the V3 position angle - V3I_YANG from north in degrees if VPARITY = -1 - Counterclockwise by the V3 position angle - V3I_YANG from north in degrees if VPARITY = 1 + Rotate hole center coordinates. + + Rotation of coordinates is: + Clockwise by the ROLL_REF + V3I_YANG from north in degrees if VPARITY = -1 + Counterclockwise by the ROLL_REF + V3I_YANG from north in degrees if VPARITY = 1 Hole center coords are in the V2, V3 plane in meters. + Notes + ----- + Nov. 2024 email discussion with Tony Sohn, Paul Goudfrooij confirmed V2/V3 coordinate + rotation back to "North up" equatorial orientation should use ROLL_REF + V3I_YANG + (= PA_APER). + Returns ------- ctrs_rot: array @@ -347,13 +356,13 @@ def mast2sky(self): # NOT used for the fringe fitting itself mask_ctrs = utils.rotate2dccw(mask_ctrs, np.pi / 2.0) vpar = self.vparity # Relative sense of rotation between Ideal xy and V2V3 - rot_ang = self.pav3 - self.v3iyang # subject to change! + rot_ang = self.roll_ref + self.v3iyang - if self.pav3 == 0.0: + if self.roll_ref == 0.0: return mask_ctrs else: # Using rotate2sccw, which rotates **vectors** CCW in a fixed coordinate system, - # so to rotate coord system CW instead of the vector, reverse sign of rotation angle. Double-check comment + # so to rotate coord system CW instead of the vector, reverse sign of rotation angle. if vpar == -1: # rotate clockwise ctrs_rot = utils.rotate2dccw(mask_ctrs, np.deg2rad(-rot_ang)) diff --git a/jwst/ami/leastsqnrm.py b/jwst/ami/leastsqnrm.py index d1c1a89524..79983d34e7 100644 --- a/jwst/ami/leastsqnrm.py +++ b/jwst/ami/leastsqnrm.py @@ -746,8 +746,8 @@ def populate_antisymmphasearray(deltaps, n=7): deltaps: 1D float array pistons between each pair of holes - n: integer - number of holes + n: integer, optional + number of holes (default=7) Returns ------- @@ -778,8 +778,8 @@ def populate_symmamparray(amps, n=7): amps: 1D float array fringe visibility between each pair of holes - n: integer - number of holes + n: integer, optional + number of holes (default=7) Returns ------- @@ -800,6 +800,43 @@ def populate_symmamparray(amps, n=7): return arr +def t3_amplitudes(amps, n=7): + """ + Populate the triple-product amplitude array + (NOT closure amplitudes) + + Parameters + ---------- + amps: 1D float array + fringe visibility between each pair of holes + + n: integer, optional + number of holes (default=7) + + Returns + ------- + cpamps: 1D float array + triple product amplitude array + """ + + arr = populate_symmamparray(amps, n=n) + + cpamps = np.zeros(int(comb(n, 3))) + + nn = 0 + for kk in range(n - 2): + for ii in range(n - kk - 2): + for jj in range(n - kk - ii - 2): + cpamps[nn + jj] = ( + arr[kk, ii + kk + 1] + * arr[ii + kk + 1, jj + ii + kk + 2] + * arr[jj + ii + kk + 2, kk] + ) + + nn += jj + 1 + + return cpamps + def redundant_cps(deltaps, n=7): """ @@ -810,8 +847,8 @@ def redundant_cps(deltaps, n=7): deltaps: 1D float array pistons between each pair of holes - n: integer - number of holes + n: integer, optional + number of holes (default=7) Returns ------- @@ -906,8 +943,8 @@ def closure_amplitudes(amps, n=7): amps: 1D float array fringe amplitudes - N: integer - number of holes + n: integer, optional + number of holes (default=7) Returns ------- @@ -934,3 +971,38 @@ def closure_amplitudes(amps, n=7): nn = nn + ll + 1 return cas + +def q4_phases(deltaps, n=7): + """ + Calculate phases for each set of 4 holes + + Parameters + ---------- + deltaps: 1D float array + pistons between each pair of holes + + n: integer, optional + number of holes (default=7) + + Returns + ------- + quad_phases: 1D float array + quad phases + """ + arr = populate_antisymmphasearray(deltaps, n=n) # fringe phase array + nn = 0 + quad_phases = np.zeros(int(comb(n, 4))) + + for ii in range(n - 3): + for jj in range(n - ii - 3): + for kk in range(n - jj - ii - 3): + for ll in range(n - jj - ii - kk - 3): + quad_phases[nn + ll] = ( + arr[ii, jj + ii + 1] + + arr[ll + ii + jj + kk + 3, kk + jj + ii + 2] + - arr[ii, kk + ii + jj + 2] + - arr[jj + ii + 1, ll + ii + jj + kk + 3] + ) + nn = nn + ll + 1 + + return quad_phases diff --git a/jwst/ami/lg_model.py b/jwst/ami/lg_model.py index 20ee86f33a..2981baaa8a 100755 --- a/jwst/ami/lg_model.py +++ b/jwst/ami/lg_model.py @@ -396,7 +396,9 @@ def fit_image( self.fringeamp, self.fringephase = leastsqnrm.tan2visibilities(self.soln) self.fringepistons = utils.fringes2pistons(self.fringephase, len(self.ctrs)) self.redundant_cps = leastsqnrm.redundant_cps(self.fringephase, n=self.N) + self.t3_amplitudes = leastsqnrm.t3_amplitudes(self.fringeamp, n=self.N) # RC 8/24 self.redundant_cas = leastsqnrm.closure_amplitudes(self.fringeamp, n=self.N) + self.q4_phases = leastsqnrm.q4_phases(self.fringephase, n=self.N) # RC 8/24 def create_modelpsf(self): """ diff --git a/jwst/ami/oifits.py b/jwst/ami/oifits.py index e0bad2e233..b373babd24 100644 --- a/jwst/ami/oifits.py +++ b/jwst/ami/oifits.py @@ -8,16 +8,20 @@ CalibOifits class: takes two AmiOIModel datamodels and produces a final calibrated datamodel. """ import numpy as np - from scipy.special import comb from astropy.stats import sigma_clipped_stats from astropy.time.core import Time +import logging from stdatamodels.jwst import datamodels +from . import leastsqnrm + +log = logging.getLogger(__name__) +log.setLevel(logging.DEBUG) class RawOifits: - def __init__(self, fringefitter, method="median"): + def __init__(self, fringefitter, method="mean"): """ Class to store AMI data in the format required to write out to OIFITS files Angular quantities of input are in radians from fringe fitting; converted to degrees for saving. @@ -28,7 +32,7 @@ def __init__(self, fringefitter, method="median"): Object containing nrm_list attribute (list of nrm objects) and other info needed for OIFITS files method: string - Method to average observables: mean or median. Default median. + Method to average observables: mean or median. Default mean. Notes ----- @@ -41,44 +45,262 @@ def __init__(self, fringefitter, method="median"): self.nslices = len(self.fringe_fitter.nrm_list) # n ints self.n_baselines = int(comb(self.n_holes, 2)) # 21 self.n_closure_phases = int(comb(self.n_holes, 3)) # 35 - self.n_closure_amplitudes = int(comb(self.n_holes, 4)) + self.n_closure_amplitudes = int(comb(self.n_holes, 4)) # also 35 self.method = method + if self.method not in ["mean", "median", "multi"]: + log.warning('"method" for saving OIFITS file must be one of ["mean","median","multi"]. Defaulting to use mean!') + self.method = "mean" self.ctrs_eqt = self.fringe_fitter.instrument_data.ctrs_eqt self.ctrs_inst = self.fringe_fitter.instrument_data.ctrs_inst - self.pa = ( - self.fringe_fitter.instrument_data.pav3 - ) # header pav3, not including v3i_yang?? self.bholes, self.bls = self._makebaselines() self.tholes, self.tuv = self._maketriples_all() + self.qholes, self.quads = self._makequads_all() def make_obsarrays(self): """ Make arrays of observables of the correct shape for saving to datamodels """ - # arrays of observables, (nslices,nobservables) shape. + # empty arrays of observables, (nslices,nobservables) shape. self.fringe_phases = np.zeros((self.nslices, self.n_baselines)) self.fringe_amplitudes = np.zeros((self.nslices, self.n_baselines)) self.closure_phases = np.zeros((self.nslices, self.n_closure_phases)) + self.t3_amplitudes = np.zeros((self.nslices, self.n_closure_phases)) + self.q4_phases = np.zeros((self.nslices, self.n_closure_amplitudes)) self.closure_amplitudes = np.zeros((self.nslices, self.n_closure_amplitudes)) self.pistons = np.zeros((self.nslices, self.n_holes)) # model parameters self.solns = np.zeros((self.nslices, 44)) + # populate with each integration's observables for i, nrmslc in enumerate(self.fringe_fitter.nrm_list): - self.fringe_phases[i, :] = np.rad2deg(nrmslc.fringephase) # FPs in degrees + self.fringe_phases[i, :] = nrmslc.fringephase # FPs in radians self.fringe_amplitudes[i, :] = nrmslc.fringeamp - self.closure_phases[i, :] = np.rad2deg(nrmslc.redundant_cps) # CPs in degrees + self.closure_phases[i, :] = nrmslc.redundant_cps # CPs in radians + self.t3_amplitudes[i, :] = nrmslc.t3_amplitudes + self.q4_phases[i, :] = nrmslc.q4_phases # quad phases in radians self.closure_amplitudes[i, :] = nrmslc.redundant_cas - self.pistons[i, :] = np.rad2deg( - nrmslc.fringepistons - ) # segment pistons in degrees + self.pistons[i, :] = nrmslc.fringepistons # segment pistons in radians self.solns[i, :] = nrmslc.soln self.fringe_amplitudes_squared = self.fringe_amplitudes ** 2 # squared visibilities + def rotate_matrix(self, cov_mat, theta): + """ + Rotate a covariance matrix by an angle. + + Parameters + ---------- + cov_mat: array + The matrix to be rotated + theta: float + Angle by which to rotate the matrix (radians) + + Returns + ------- + cv_rotated: array + The rotated matrix + + """ + c, s = np.cos(theta), np.sin(theta) + r_mat = [[c, -s], + [s, c]] + # coordinate rotation from real/imaginary to absolute value/phase (modulus/argument) + return np.linalg.multi_dot([np.transpose(r_mat), cov_mat, r_mat]) + + + def average_observables(self, averfunc): + """ + Average all the observables. + + Calculate covariance matrices between fringe amplitudes/fringe phases, + and between triple product amps/closure phases, and closure amplitudes/quad phases. + Convert r, theta (modulus, phase) to x,y. Calculate cov(x,y). Rotate resulting + 2x2 matrix back to r, theta. Take sqrt of relevant covariance element to be error. + This must be done with phases in radians. + + Parameters + ---------- + averfunc: function + np.mean (default) or np.median + + Returns + ------- + avg_sqv: array + Averaged squared visibilites + err_sqv: array + Standard error of the mean of averaged squared visibilities + avg_fa: array + Averaged fringe (visibility) amplitudes + err_fa: array + Standard error of the mean of averaged fringe (visibility) amplitudes + avg_fp: array + Averaged fringe phases (rad) + err_fp: array + Standard error of the mean of averaged fringe phases (rad) + avg_cp: array + Averaged closure phases (rad) + err_cp: array + Standard error of the mean of averaged closure phases (rad) + avg_t3amp: array + Averaged triple amplitudes + err_t3amp: array + Standard error of the mean of averaged triple amplitudes + avg_ca: array + Averaged closure amplitudes + err_ca: array + Standard error of the mean of averaged closure amplitudes + avg_q4phi: array + Averaged quad phases + err_q4phi: array + Standard error of the mean of averaged quad phases + avg_pist: array + Averaged segment pistons + err_pist: array + Standard error of the mean of averaged segment pistons + + """ + covmats_fringes, covmats_triples, covmats_quads = self.observable_covariances(averfunc) + + if self.method == 'mean': + avg_fa, _, std_fa = sigma_clipped_stats(self.fringe_amplitudes, axis=0) + avg_fp, _, std_fp = sigma_clipped_stats(self.fringe_phases, axis=0) + avg_sqv, _, std_sqv = sigma_clipped_stats(self.fringe_amplitudes**2, axis=0) + avg_pist, _, err_pist = sigma_clipped_stats(self.pistons, axis=0) + else: # median + _, avg_fa, std_fa = sigma_clipped_stats(self.fringe_amplitudes, axis=0) # 21. std_fa is just for comparing to covariance + _, avg_fp, std_fp = sigma_clipped_stats(self.fringe_phases, axis=0) # 21 + _, avg_sqv, std_sqv = sigma_clipped_stats(self.fringe_amplitudes**2, axis=0) + _, avg_pist, err_pist = sigma_clipped_stats(self.pistons, axis=0) + + err_pist = err_pist/np.sqrt(self.nslices) # standard error of the mean + err_fa, err_fp = self.err_from_covmat(covmats_fringes) + + # calculate squared visibility (fringe) amplitude uncertainties correctly + err_sqv = (2 * avg_fa * err_fa) / np.sqrt(self.nslices) + + # calculate triple and quad quantities from **averaged** fringe amps and phases + avg_t3amp = leastsqnrm.t3_amplitudes(avg_fa, n=self.n_holes) + avg_cp = leastsqnrm.redundant_cps(avg_fp, n=self.n_holes) + err_t3amp, err_cp = self.err_from_covmat(covmats_triples) + + avg_ca = leastsqnrm.closure_amplitudes(avg_fa, n=self.n_holes) + avg_q4phi = leastsqnrm.q4_phases(avg_fp, n=self.n_holes) + err_ca, err_q4phi = self.err_from_covmat(covmats_quads) + + return (avg_sqv, err_sqv, avg_fa, err_fa, + avg_fp, err_fp, avg_cp, err_cp, + avg_t3amp, err_t3amp, avg_ca, err_ca, + avg_q4phi, err_q4phi, avg_pist, err_pist) + + def err_from_covmat(self, covmatlist): + """ + Derive observable errors from their covariance matrices. + + Return sqrt of [0,0] and [1,1] elements of each of a list of covariance matrices, + divided by sqrt(N_ints), for use as observable errors (standard error of the mean). + If using median, error calculation is questionable because this is NOT the standard + error of the median. + + Parameters + ---------- + covmatlist: array + array of covariance matrices for each baseline/triple/quad + shape e.g. (21,2,2) or (35,2,2) + + Returns + ------- + err_00: array + standard errors of the mean of the first observable. shape e.g. (21) + err_11: array + standard errors of the mean of the second observable. shape e.g. (21) + """ + err_00 = np.sqrt(np.array([covmat[0,0] for covmat in covmatlist]))/np.sqrt(self.nslices) + err_11 = np.sqrt(np.array([covmat[1,1] for covmat in covmatlist]))/np.sqrt(self.nslices) + + return err_00, err_11 + + def observable_covariances(self, averfunc): + """ + Calculate covariance matrices from each pair of observables. + + For each baseline/triple/quad, calculate covariance between each + fringe amplitude/phase quantity. + + Parameters + ---------- + averfunc: function + np.mean (default) or np.median + + Returns + ------- + cov_mat_fringes: array + Array of 21 covariance matrices for fringes (amplitudes, phases) + cov_mat_triples: array + Array of 35 covariance matrices for triples (t3 amplitudes, closure phases) + cov_mat_quads: array + Array of 35 covariance matrices for quads (closure amplitudes, quad phases) + + """ + # loop over 21 baselines + cov_mat_fringes = [] + # These operate on all slices at once, not on already-averaged integrations + for bl in np.arange(self.n_baselines): + fringeamps = self.fringe_amplitudes[:,bl] + fringephases = self.fringe_phases[:,bl] + covmat = self.cov_r_theta(fringeamps, fringephases, averfunc) + cov_mat_fringes.append(covmat) + # loop over 35 triples + cov_mat_triples = [] + for triple in np.arange(self.n_closure_phases): + tripamp = self.t3_amplitudes[:,triple] + triphase = self.closure_phases[:,triple] + covmat = self.cov_r_theta(tripamp, triphase, averfunc) + cov_mat_triples.append(covmat) + # loop over 35 quads + cov_mat_quads = [] + for quad in np.arange(self.n_closure_amplitudes): + quadamp = self.closure_amplitudes[:,quad] + quadphase = self.q4_phases[:,quad] + covmat = self.cov_r_theta(quadamp, quadphase, averfunc) + cov_mat_quads.append(covmat) + + # covmats to be written to oifits. store in rawoifits object? TBD + # lists of cov mats have shape e.g. (21, 2, 2) or (35, 2, 2) + + return np.array(cov_mat_fringes), np.array(cov_mat_triples), np.array(cov_mat_quads) + + + def cov_r_theta(self, rr, theta, averfunc): + """ + Calculate covariance in polar coordinates. + + Calculate covariance in x, y coordinates, then rotate covariance matrix + by **average** phase (over integrations) to get matrix in (r,theta). + + Parameters + ---------- + rr: array + complex number modulus + theta: array + complex number phase + averfunc: function + np.mean (default) or np.median + + Returns + ------- + cov_mat_r_theta: array (2,2) + Covariance matrix in r, theta coordinates + + """ + xx = rr * np.cos(theta) + yy = rr * np.sin(theta) + cov_mat_xy = np.cov(xx, yy) + return self.rotate_matrix(cov_mat_xy, averfunc(theta)) + + def make_oifits(self): """ Perform final manipulations of observable arrays, calculate uncertainties, and @@ -114,16 +336,14 @@ def make_oifits(self): flagVis = [False] * self.n_baselines flagT3 = [False] * self.n_closure_phases - # do the things done by populate_nrm here - # average or don't, and get uncertainties + # Average observables (or don't), and get uncertainties # Unwrap phases shift2pi = np.zeros(self.closure_phases.shape) shift2pi[self.closure_phases >= 6] = 2 * np.pi shift2pi[self.closure_phases <= -6] = -2 * np.pi self.closure_phases -= shift2pi - if self.method not in ["mean", "median", "multi"]: - self.method = "median" + # Now we are setting up the observables to be written out to OIFITS # set these as attributes (some may exist and be overwritten) if self.method == "multi": self.vis2 = self.fringe_amplitudes_squared.T @@ -134,35 +354,43 @@ def make_oifits(self): self.e_visphi = np.zeros(self.visphi.shape) self.closure_phases = self.closure_phases.T self.e_cp = np.zeros(self.closure_phases.shape) + self.t3amp = self.t3_amplitudes.T + self.e_t3amp = np.zeros(self.t3amp.shape) + self.q4phi = self.q4_phases.T + self.e_q4phi = np.zeros(self.q4_phases.shape) self.camp = self.closure_amplitudes.T self.e_camp = np.zeros(self.camp.shape) self.pist = self.pistons.T self.e_pist = np.zeros(self.pist.shape) - # apply sigma-clipping to uncertainties - # sigma_clipped_stats returns mean, median, stddev. nsigma=3, niters=5 - elif self.method == "median": - _, self.vis2, self.e_vis2 = sigma_clipped_stats(self.fringe_amplitudes_squared, axis=0) - _, self.visamp, self.e_visamp = sigma_clipped_stats(self.fringe_amplitudes, axis=0) - _, self.visphi, self.e_visphi = sigma_clipped_stats(self.fringe_phases, axis=0) - _, self.closure_phases, self.e_cp = sigma_clipped_stats(self.closure_phases, axis=0) - _, self.camp, self.e_camp = sigma_clipped_stats(self.closure_amplitudes, axis=0) - _, self.pist, self.e_pist = sigma_clipped_stats(self.pistons, axis=0) - - else: # take the mean - self.vis2, _, self.e_vis2 = sigma_clipped_stats(self.fringe_amplitudes_squared, axis=0) - self.visamp, _, self.e_visamp = sigma_clipped_stats(self.fringe_amplitudes, axis=0) - self.visphi, _, self.e_visphi = sigma_clipped_stats(self.fringe_phases, axis=0) - self.closure_phases, _, self.e_cp = sigma_clipped_stats(self.closure_phases, axis=0) - self.camp, _, self.e_camp = sigma_clipped_stats(self.closure_amplitudes, axis=0) - self.pist, _, self.e_pist = sigma_clipped_stats(self.pistons, axis=0) + + elif self.method == "mean": + (self.vis2, self.e_vis2, self.visamp, self.e_visamp, + self.visphi, self.e_visphi, self.closure_phases, self.e_cp, + self.t3amp, self.e_t3amp, self.camp, self.e_camp, + self.q4phi, self.e_q4phi, self.pist, self.e_pist) = self.average_observables(np.mean) + + else: # take the median + (self.vis2, self.e_vis2, self.visamp, self.e_visamp, + self.visphi, self.e_visphi, self.closure_phases, self.e_cp, + self.t3amp, self.e_t3amp, self.camp, self.e_camp, + self.q4phi, self.e_q4phi, self.pist, self.e_pist) = self.average_observables(np.median) + + # Convert angular quantities from radians to degrees + self.visphi = np.rad2deg(self.visphi) + self.e_visphi = np.rad2deg(self.e_visphi) + self.closure_phases = np.rad2deg(self.closure_phases) + self.e_cp = np.rad2deg(self.e_cp) + self.q4phi = np.rad2deg(self.q4phi) + self.e_q4phi = np.rad2deg(self.e_q4phi) + self.pist = np.rad2deg(self.pist) + self.e_pist = np.rad2deg(self.e_pist) # prepare arrays for OI_ARRAY ext self.staxy = instrument_data.ctrs_inst - N_ap = len(self.staxy) - tel_name = ["A%i" % x for x in np.arange(N_ap) + 1] + tel_name = ["A%i" % x for x in np.arange(self.n_holes) + 1] sta_name = tel_name - diameter = [0] * N_ap + diameter = [0] * self.n_holes staxyz = [] for x in self.staxy: @@ -170,105 +398,108 @@ def make_oifits(self): line = [a[0], a[1], 0] staxyz.append(line) - sta_index = np.arange(N_ap) + 1 + sta_index = np.arange(self.n_holes) + 1 pscale = instrument_data.pscale_mas / 1000.0 # arcsec - isz = self.fringe_fitter.scidata.shape[ - 1 - ] # Size of the image to extract NRM data - fov = [pscale * isz] * N_ap - fovtype = ["RADIUS"] * N_ap + # Size of the image to extract NRM data + isz = self.fringe_fitter.scidata.shape[1] + fov = [pscale * isz] * self.n_holes + fovtype = ["RADIUS"] * self.n_holes - m = datamodels.AmiOIModel() - self.init_oimodel_arrays(m) + oim = datamodels.AmiOIModel() + self.init_oimodel_arrays(oim) # primary header keywords - m.meta.telescope = instrument_data.telname - m.meta.origin = "STScI" - m.meta.instrument.name = instrument_data.instrument - m.meta.program.pi_name = instrument_data.pi_name - m.meta.target.proposer_name = instrument_data.proposer_name - m.meta.observation.date = observation_date.fits - m.meta.oifits.array_name = instrument_data.arrname - m.meta.oifits.instrument_mode = instrument_data.pupil + oim.meta.telescope = instrument_data.telname + oim.meta.origin = "STScI" + oim.meta.instrument.name = instrument_data.instrument + oim.meta.program.pi_name = instrument_data.pi_name + oim.meta.target.proposer_name = instrument_data.proposer_name + oim.meta.observation.date = observation_date.fits + oim.meta.oifits.array_name = instrument_data.arrname + oim.meta.oifits.instrument_mode = instrument_data.pupil + + oim.meta.ami.roll_ref = instrument_data.roll_ref + oim.meta.ami.v3yangle = instrument_data.v3iyang + oim.meta.ami.vparity = instrument_data.vparity # oi_array extension data - m.array["TEL_NAME"] = tel_name - m.array["STA_NAME"] = sta_name - m.array["STA_INDEX"] = sta_index - m.array["DIAMETER"] = diameter - m.array["STAXYZ"] = staxyz - m.array["FOV"] = fov - m.array["FOVTYPE"] = fovtype - m.array["CTRS_EQT"] = instrument_data.ctrs_eqt - m.array["PISTONS"] = self.pist - m.array["PIST_ERR"] = self.e_pist + oim.array["TEL_NAME"] = tel_name + oim.array["STA_NAME"] = sta_name + oim.array["STA_INDEX"] = sta_index + oim.array["DIAMETER"] = diameter + oim.array["STAXYZ"] = staxyz + oim.array["FOV"] = fov + oim.array["FOVTYPE"] = fovtype + oim.array["CTRS_EQT"] = instrument_data.ctrs_eqt + oim.array["PISTONS"] = self.pist + oim.array["PIST_ERR"] = self.e_pist # oi_target extension data - m.target['TARGET_ID'] = [1] - m.target['TARGET'] = instrument_data.objname - m.target['RAEP0'] = instrument_data.ra - m.target['DECEP0'] = instrument_data.dec - m.target['EQUINOX'] = [2000] - m.target['RA_ERR'] = instrument_data.ra_uncertainty - m.target['DEC_ERR'] = instrument_data.dec_uncertainty - m.target['SYSVEL'] = [0] - m.target['VELTYP'] = ['UNKNOWN'] - m.target['VELDEF'] = ['OPTICAL'] - m.target['PMRA'] = instrument_data.pmra - m.target['PMDEC'] = instrument_data.pmdec - m.target['PMRA_ERR'] = [0] - m.target['PMDEC_ERR'] = [0] - m.target['PARALLAX'] = [0] - m.target['PARA_ERR'] = [0] - m.target['SPECTYP'] = ['UNKNOWN'] + oim.target['TARGET_ID'] = [1] + oim.target['TARGET'] = instrument_data.objname + oim.target['RAEP0'] = instrument_data.ra + oim.target['DECEP0'] = instrument_data.dec + oim.target['EQUINOX'] = [2000] + oim.target['RA_ERR'] = instrument_data.ra_uncertainty + oim.target['DEC_ERR'] = instrument_data.dec_uncertainty + oim.target['SYSVEL'] = [0] + oim.target['VELTYP'] = ['UNKNOWN'] + oim.target['VELDEF'] = ['OPTICAL'] + oim.target['PMRA'] = instrument_data.pmra + oim.target['PMDEC'] = instrument_data.pmdec + oim.target['PMRA_ERR'] = [0] + oim.target['PMDEC_ERR'] = [0] + oim.target['PARALLAX'] = [0] + oim.target['PARA_ERR'] = [0] + oim.target['SPECTYP'] = ['UNKNOWN'] # oi_vis extension data - m.vis['TARGET_ID'] = 1 - m.vis['TIME'] = 0 - m.vis['MJD'] = observation_date.mjd - m.vis['INT_TIME'] = instrument_data.itime - m.vis['VISAMP'] = self.visamp - m.vis['VISAMPERR'] = self.e_visamp - m.vis['VISPHI'] = self.visphi - m.vis['VISPHIERR'] = self.e_visphi - m.vis['UCOORD'] = ucoord - m.vis['VCOORD'] = vcoord - m.vis['STA_INDEX'] = self._format_STAINDEX_V2(self.bholes) - m.vis['FLAG'] = flagVis + oim.vis['TARGET_ID'] = 1 + oim.vis['TIME'] = 0 + oim.vis['MJD'] = observation_date.mjd + oim.vis['INT_TIME'] = instrument_data.itime + oim.vis['VISAMP'] = self.visamp + oim.vis['VISAMPERR'] = self.e_visamp + oim.vis['VISPHI'] = self.visphi + oim.vis['VISPHIERR'] = self.e_visphi + oim.vis['UCOORD'] = ucoord + oim.vis['VCOORD'] = vcoord + oim.vis['STA_INDEX'] = self._format_STAINDEX_V2(self.bholes) + oim.vis['FLAG'] = flagVis # oi_vis2 extension data - m.vis2['TARGET_ID'] = 1 - m.vis2['TIME'] = 0 - m.vis2['MJD'] = observation_date.mjd - m.vis2['INT_TIME'] = instrument_data.itime - m.vis2['VIS2DATA'] = self.vis2 - m.vis2['VIS2ERR'] = self.e_vis2 - m.vis2['UCOORD'] = ucoord - m.vis2['VCOORD'] = vcoord - m.vis2['STA_INDEX'] = self._format_STAINDEX_V2(self.bholes) - m.vis2['FLAG'] = flagVis + oim.vis2['TARGET_ID'] = 1 + oim.vis2['TIME'] = 0 + oim.vis2['MJD'] = observation_date.mjd + oim.vis2['INT_TIME'] = instrument_data.itime + oim.vis2['VIS2DATA'] = self.vis2 + oim.vis2['VIS2ERR'] = self.e_vis2 + oim.vis2['UCOORD'] = ucoord + oim.vis2['VCOORD'] = vcoord + oim.vis2['STA_INDEX'] = self._format_STAINDEX_V2(self.bholes) + oim.vis2['FLAG'] = flagVis # oi_t3 extension data - m.t3['TARGET_ID'] = 1 - m.t3['TIME'] = 0 - m.t3['MJD'] = observation_date.mjd - m.t3['T3AMP'] = self.camp - m.t3['T3AMPERR'] = self.e_camp - m.t3['T3PHI'] = self.closure_phases - m.t3['T3PHIERR'] = self.e_cp - m.t3['U1COORD'] = u1coord - m.t3['V1COORD'] = v1coord - m.t3['U2COORD'] = u2coord - m.t3['V2COORD'] = v2coord - m.t3['STA_INDEX'] = self._format_STAINDEX_T3(self.tholes) - m.t3['FLAG'] = flagT3 + oim.t3['TARGET_ID'] = 1 + oim.t3['TIME'] = 0 + oim.t3['MJD'] = observation_date.mjd + oim.t3['T3AMP'] = self.t3amp + oim.t3['T3AMPERR'] = self.e_t3amp + oim.t3['T3PHI'] = self.closure_phases + oim.t3['T3PHIERR'] = self.e_cp + oim.t3['U1COORD'] = u1coord + oim.t3['V1COORD'] = v1coord + oim.t3['U2COORD'] = u2coord + oim.t3['V2COORD'] = v2coord + oim.t3['STA_INDEX'] = self._format_STAINDEX_T3(self.tholes) + oim.t3['FLAG'] = flagT3 # oi_wavelength extension data - m.wavelength["EFF_WAVE"] = wl - m.wavelength["EFF_BAND"] = e_wl + oim.wavelength["EFF_WAVE"] = wl + oim.wavelength["EFF_BAND"] = e_wl - return m + return oim def init_oimodel_arrays(self, oimodel): """ @@ -381,29 +612,22 @@ def _maketriples_all(self): ------- tarray: integer array Triple hole indices (0-indexed), - float array of two uv vectors in all triangles + float array of two uv vectors in all triangles """ tlist = [] + uvlist = [] for i in range(self.n_holes): for j in range(self.n_holes): for k in range(self.n_holes): if i < j and j < k: tlist.append((i, j, k)) + uvlist.append( + ( + self.ctrs_eqt[i] - self.ctrs_eqt[j], + self.ctrs_eqt[j] - self.ctrs_eqt[k], + ) + ) tarray = np.array(tlist).astype(int) - - tname = [] - uvlist = [] - # foreach row of 3 elts... - for triple in tarray: - tname.append("{0:d}_{1:d}_{2:d}".format(triple[0], triple[1], triple[2])) - - uvlist.append( - ( - self.ctrs_eqt[triple[0]] - self.ctrs_eqt[triple[1]], - self.ctrs_eqt[triple[1]] - self.ctrs_eqt[triple[2]], - ) - ) - return tarray, np.array(uvlist) def _makebaselines(self): @@ -425,6 +649,29 @@ def _makebaselines(self): bllist.append(self.ctrs_eqt[i] - self.ctrs_eqt[j]) return np.array(blist).astype(int), np.array(bllist) + def _makequads_all(self): + """ + Calculate all four-hole combinations (quads) + + Returns + ------- + qarray: int array of four-hole quads (0-based) + uvwlist: numpy array of u, v, w vectors for each quad + """ + qlist = [] + uvwlist = [] + for i in range(self.n_holes): + for j in range(self.n_holes): + for k in range(self.n_holes): + for q in range(self.n_holes): + if i < j and j < k and k < q: + qlist.append((i, j, k, q)) + uvwlist.append((self.ctrs_eqt[i] - self.ctrs_eqt[j], + self.ctrs_eqt[j] - self.ctrs_eqt[k], + self.ctrs_eqt[k] - self.ctrs_eqt[q])) + qarray = np.array(qlist).astype(int) + return qarray, np.array(uvwlist) + def _format_STAINDEX_T3(self, tab): """ Converts sta_index to save oifits T3 in the appropriate format @@ -581,6 +828,6 @@ def calibrate(self): # add calibrated header keywords calname = self.caloimodel.meta.target.proposer_name # name of calibrator star - self.calib_oimodel.meta.oifits.calib = calname + self.calib_oimodel.meta.ami.calibrator_object_id = calname return self.calib_oimodel diff --git a/jwst/ami/tests/test_ami_interface.py b/jwst/ami/tests/test_ami_interface.py index 2da7d16eb3..23e765a343 100644 --- a/jwst/ami/tests/test_ami_interface.py +++ b/jwst/ami/tests/test_ami_interface.py @@ -48,7 +48,7 @@ def example_model(mock_nrm_reference_file): model.meta.program.pi_name = "someone" model.meta.target.catalog_name = "" model.meta.visit.start_time = "2022-06-05 12:15:41.5020000" - model.meta.pointing.pa_v3 = 171.8779402866089 + model.meta.wcsinfo.roll_ref = 171.8779402866089 model.meta.wcsinfo.v3yangle = 0.56126717 model.meta.filename = "test_calints.fits" model.meta.instrument.pupil = "NRM" @@ -63,7 +63,6 @@ def test_ami_analyze_even_oversample_fail(example_model, oversample): AmiAnalyzeStep.call(example_model, oversample=oversample) -@pytest.mark.skip("reference files are not currently used") def test_ami_analyze_no_reffile_fail(monkeypatch, example_model): """Make sure that ami_analyze fails if no throughput reffile is available""" @@ -71,7 +70,7 @@ def mockreturn(input_model, reftype, observatory=None, asn_exptypes=None): return "N/A" monkeypatch.setattr(stpipe.crds_client, 'get_reference_file', mockreturn) - with pytest.raises(RuntimeError, match="No throughput reference file found."): + with pytest.raises(RuntimeError, match="No THROUGHPUT reference file found."): AmiAnalyzeStep.call(example_model) diff --git a/jwst/ami/utils.py b/jwst/ami/utils.py index 700f87a2c9..a9fa5e1712 100644 --- a/jwst/ami/utils.py +++ b/jwst/ami/utils.py @@ -1322,16 +1322,12 @@ def handle_bandpass(bandpass, throughput_model): Array of weights, wavelengths used to generate model """ if bandpass is not None: - log.info( - "User-defined bandpass provided: OVERWRITING ALL NIRISS-SPECIFIC FILTER/BANDPASS VARIABLES" - ) # bandpass can be user-defined synphot object or appropriate array if isinstance(bandpass, synphot.spectrum.SpectralElement): log.info("User-defined synphot spectrum provided") wl, wt = bandpass._get_arrays(bandpass.waveset) bandpass = np.array((wt, wl)).T else: - log.info("User-defined bandpass array provided") bandpass = np.array(bandpass) else: diff --git a/jwst/assign_wcs/nirspec.py b/jwst/assign_wcs/nirspec.py index 7902d6c3bc..26149cc262 100644 --- a/jwst/assign_wcs/nirspec.py +++ b/jwst/assign_wcs/nirspec.py @@ -726,8 +726,7 @@ def get_open_msa_slits(prog_id, msa_file, msa_metadata_id, dither_position, quadrant = slitlet_rows[0]['shutter_quadrant'] ycen = j xcen = slitlet_rows[0]['shutter_row'] # grab the first as they are all the same - shutter_id = xcen + (ycen - 1) * 365 # shutter numbers in MSA file are 1-indexed - + shutter_id = np.int64(xcen) + (np.int64(ycen) - 1) * 365 # shutter numbers in MSA file are 1-indexed # Background slits all have source_id=0 in the msa_file, # so assign a unique id based on the slitlet_id source_id = slitlet_id @@ -750,7 +749,7 @@ def get_open_msa_slits(prog_id, msa_file, msa_metadata_id, dither_position, np.nan_to_num(s['estimated_source_in_shutter_x'], nan=0.5), np.nan_to_num(s['estimated_source_in_shutter_y'], nan=0.5)) for s in slitlet_rows if s['background'] == 'N'][0] - shutter_id = xcen + (ycen - 1) * 365 # shutter numbers in MSA file are 1-indexed + shutter_id = np.int64(xcen) + (np.int64(ycen) - 1) * 365 # shutter numbers in MSA file are 1-indexed # y-size jmin = min([s['shutter_column'] for s in slitlet_rows]) diff --git a/jwst/jump/jump.py b/jwst/jump/jump.py index 92f70de2fb..febae44348 100644 --- a/jwst/jump/jump.py +++ b/jwst/jump/jump.py @@ -28,7 +28,8 @@ def run_detect_jumps(output_model, gain_model, readnoise_model, minimum_sigclip_groups=100, only_use_ints=True, mask_snowball_persist_next_int=True, - snowball_time_masked_next_int=250 + snowball_time_masked_next_int=250, + max_shower_amplitude=4 ): # Runs `detect_jumps` in stcal @@ -43,6 +44,8 @@ def run_detect_jumps(output_model, gain_model, readnoise_model, after_jump_flag_n2 = int(after_jump_flag_time2 // gtime) grps_masked_after_shower = int(time_masked_after_shower // gtime) snowball_grps_masked_next_int = int(snowball_time_masked_next_int // gtime) + # Likewise, convert a max MIRI shower amplitude in DN/s to DN/group + max_shower_amplitude = max_shower_amplitude * gtime # Get 2D gain and read noise values from their respective models if reffile_utils.ref_matches_sci(output_model, gain_model): gain_2d = gain_model.data @@ -85,7 +88,8 @@ def run_detect_jumps(output_model, gain_model, readnoise_model, minimum_sigclip_groups=minimum_sigclip_groups, only_use_ints=only_use_ints, mask_persist_grps_next_int = mask_snowball_persist_next_int, - persist_grps_flagged = snowball_grps_masked_next_int + persist_grps_flagged = snowball_grps_masked_next_int, + max_shower_amplitude = max_shower_amplitude ) diff --git a/jwst/jump/jump_step.py b/jwst/jump/jump_step.py index 86f2b33f23..7eec5c1780 100755 --- a/jwst/jump/jump_step.py +++ b/jwst/jump/jump_step.py @@ -38,6 +38,7 @@ class JumpStep(Step): mask_snowball_core_next_int = boolean(default=True) # Flag saturated cores of snowballs in the next integration? snowball_time_masked_next_int = integer(default=4000) # Time in seconds over which saturated cores are flagged in next integration find_showers = boolean(default=False) # Apply MIRI shower flagging? + max_shower_amplitude = float(default=4) # Maximum MIRI shower amplitude in DN/s extend_snr_threshold = float(default=1.2) # The SNR minimum for detection of extended showers in MIRI extend_min_area = integer(default=90) # Min area of emission after convolution for the detection of showers extend_inner_radius = float(default=1) # Inner radius of the ring_2D_kernel used for convolution @@ -119,6 +120,7 @@ def process(self, step_input): min_sat_radius_extend=self.min_sat_radius_extend, sat_required_snowball=sat_required_snowball, sat_expand=self.sat_expand * 2, expand_large_events=expand_large_events, find_showers=self.find_showers, + max_shower_amplitude=self.max_shower_amplitude, edge_size=self.edge_size, extend_snr_threshold=self.extend_snr_threshold, extend_min_area=self.extend_min_area, extend_inner_radius=self.extend_inner_radius, diff --git a/jwst/lib/tests/test_engdb_mast.py b/jwst/lib/tests/test_engdb_mast.py index faa6d242b8..66a17df001 100644 --- a/jwst/lib/tests/test_engdb_mast.py +++ b/jwst/lib/tests/test_engdb_mast.py @@ -67,22 +67,23 @@ def test_get_records(engdb): [ ({}, [-0.7914494276, -0.7914494276, -0.7914494276, -0.791449368]), ({'include_obstime': True}, - [EngDB_Value(obstime=Time('2022-02-02T22:24:58.053', scale='utc', format='isot'), value=-0.7914494276), - EngDB_Value(obstime=Time('2022-02-02T22:24:58.309', scale='utc', format='isot'), value=-0.7914494276), - EngDB_Value(obstime=Time('2022-02-02T22:24:58.565', scale='utc', format='isot'), value=-0.7914494276), - EngDB_Value(obstime=Time('2022-02-02T22:24:58.821', scale='utc', format='isot'), value=-0.791449368)]), + [EngDB_Value(obstime=Time(59612.9340052431, format='mjd'), value=-0.7914494276), + EngDB_Value(obstime=Time(59612.9340082060, format='mjd'), value=-0.7914494276), + EngDB_Value(obstime=Time(59612.9340111690, format='mjd'), value=-0.7914494276), + EngDB_Value(obstime=Time(59612.9340141319, format='mjd'), value=-0.791449368)]), ({'include_obstime': True, 'zip_results': False}, EngDB_Value( - obstime=[Time('2022-02-02T22:24:58.053', scale='utc', format='isot'), - Time('2022-02-02T22:24:58.309', scale='utc', format='isot'), - Time('2022-02-02T22:24:58.565', scale='utc', format='isot'), - Time('2022-02-02T22:24:58.821', scale='utc', format='isot')], + obstime=[ + Time(59612.9340052431, format='mjd'), + Time(59612.9340082060, format='mjd'), + Time(59612.9340111690, format='mjd'), + Time(59612.9340141319, format='mjd')], value=[-0.7914494276, -0.7914494276, -0.7914494276, -0.791449368])), ({'include_bracket_values': True}, [-0.791449368, -0.7914494276, -0.7914494276, -0.7914494276, -0.791449368, -0.791449368]) ]) def test_get_values(engdb, pars, expected): values = engdb.get_values(*QUERY, **pars) - assert str(values) == str(expected) + assert values == expected def test_negative_aliveness(): diff --git a/jwst/master_background/master_background_mos_step.py b/jwst/master_background/master_background_mos_step.py index 78a6bcfb96..d5ad86b711 100644 --- a/jwst/master_background/master_background_mos_step.py +++ b/jwst/master_background/master_background_mos_step.py @@ -9,6 +9,9 @@ from ..flatfield import flat_field_step from ..pathloss import pathloss_step from ..photom import photom_step +from ..pixel_replace import pixel_replace_step +from ..resample import resample_spec_step +from ..extract_1d import extract_1d_step from ..stpipe import Pipeline __all__ = ['MasterBackgroundMosStep'] @@ -62,6 +65,9 @@ class MasterBackgroundMosStep(Pipeline): 'pathloss': pathloss_step.PathLossStep, 'barshadow': barshadow_step.BarShadowStep, 'photom': photom_step.PhotomStep, + 'pixel_replace': pixel_replace_step.PixelReplaceStep, + 'resample_spec': resample_spec_step.ResampleSpecStep, + 'extract_1d': extract_1d_step.Extract1dStep, } # No need to prefetch. This will have been done by the parent step. @@ -170,7 +176,7 @@ def process(self, data): return result def set_pars_from_parent(self): - """Set substep parameters from the parents substeps""" + """Set substep parameters from the parents substeps when needed""" if not self.parent: return @@ -188,6 +194,21 @@ def set_pars_from_parent(self): del pars[par] getattr(self, step).update_pars(pars) + def _extend_bg_slits(self, pre_calibrated): + # Copy dedicated background slitlets to a temporary model + bkg_model = datamodels.MultiSlitModel() + bkg_model.update(pre_calibrated) + slits = [] + for slit in pre_calibrated.slits: + if nirspec_utils.is_background_msa_slit(slit): + self.log.info(f'Using background slitlet {slit.source_name}') + slits.append(slit) + if len(slits) == 0: + self.log.warning('No background slitlets found; skipping master bkg correction') + return None + bkg_model.slits.extend(slits) + return bkg_model + def _classify_slits(self, data): """Determine how many Slits are background and source types @@ -269,9 +290,16 @@ def _calc_master_background( master_background = user_background bkg_x1d_spectra = None else: - self.log.debug('Calculating 1D master background') - master_background, bkg_x1d_spectra = nirspec_utils.create_background_from_multislit( - pre_calibrated, sigma_clip=sigma_clip, median_kernel=median_kernel) + self.log.info('Creating MOS master background from background slitlets') + bkg_model = self._extend_bg_slits(pre_calibrated) + if bkg_model is not None: + bkg_model = self.pixel_replace.run(bkg_model) + bkg_model = self.resample_spec.run(bkg_model) + bkg_x1d_spectra = self.extract_1d.run(bkg_model) + master_background = nirspec_utils.create_background_from_multispec( + bkg_x1d_spectra, sigma_clip=sigma_clip, median_kernel=median_kernel) + else: + master_background = None if master_background is None: self.log.debug('No master background could be calculated. Returning None') return None, None, None diff --git a/jwst/master_background/nirspec_utils.py b/jwst/master_background/nirspec_utils.py index 862e6b82af..3420deca6a 100644 --- a/jwst/master_background/nirspec_utils.py +++ b/jwst/master_background/nirspec_utils.py @@ -3,7 +3,6 @@ from scipy.signal import medfilt -from stdatamodels.jwst import datamodels log = logging.getLogger(__name__) log.setLevel(logging.DEBUG) @@ -78,14 +77,14 @@ def map_to_science_slits(input_model, master_bkg): return output_model -def create_background_from_multislit(input_model, sigma_clip=3, median_kernel=1): +def create_background_from_multispec(bkg_model, sigma_clip=3, median_kernel=1): """Create a 1D master background spectrum from a set of calibrated background MOS slitlets in the input - MultiSlitModel. + MultiSpecModel. Parameters ---------- - input_model : `~jwst.datamodels.MultiSlitModel` + bkg_model : `~jwst.datamodels.MultiSpecModel` The input data model containing all slit instances. sigma_clip : None or float, optional Optional factor for sigma clipping outliers when combining background spectra. @@ -101,36 +100,11 @@ def create_background_from_multislit(input_model, sigma_clip=3, median_kernel=1) x1d: `jwst.datamodels.MultiSpecModel` The 1D extracted background spectra of the inputs. """ - from ..resample import resample_spec_step - from ..extract_1d import extract_1d_step from ..combine_1d.combine1d import combine_1d_spectra - log.info('Creating MOS master background from background slitlets') - - # Copy dedicated background slitlets to a temporary model - bkg_model = datamodels.MultiSlitModel() - bkg_model.update(input_model) - slits = [] - for slit in input_model.slits: - if is_background_msa_slit(slit): - log.info(f'Using background slitlet {slit.source_name}') - slits.append(slit) - - if len(slits) == 0: - log.warning('No background slitlets found; skipping master bkg correction') - return None - - bkg_model.slits.extend(slits) - - # Apply resample_spec and extract_1d to all background slitlets - log.info('Applying resampling and 1D extraction to background slits') - resamp = resample_spec_step.ResampleSpecStep.call(bkg_model) - x1d = extract_1d_step.Extract1dStep.call(resamp) - # Call combine_1d to combine the 1D background spectra log.info('Combining 1D background spectra into master background') - master_bkg = combine_1d_spectra( - x1d, exptime_key='exposure_time', sigma_clip=sigma_clip) + master_bkg = combine_1d_spectra(bkg_model, exptime_key='exposure_time', sigma_clip=sigma_clip) # If requested, apply a moving-median boxcar filter to the master background spectrum # Round down even kernel sizes because only odd kernel sizes are supported. @@ -150,10 +124,7 @@ def create_background_from_multislit(input_model, sigma_clip=3, median_kernel=1) kernel_size=[median_kernel] ) - del bkg_model - del resamp - - return master_bkg, x1d + return master_bkg def correct_nrs_ifu_bkg(input_model): diff --git a/jwst/master_background/tests/test_master_background_mos.py b/jwst/master_background/tests/test_master_background_mos.py index b683f49616..f0fc67042d 100644 --- a/jwst/master_background/tests/test_master_background_mos.py +++ b/jwst/master_background/tests/test_master_background_mos.py @@ -2,7 +2,7 @@ import pytest from astropy.io import fits from astropy.table import Table -from stdatamodels.jwst.datamodels import ImageModel +from stdatamodels.jwst.datamodels import ImageModel, MultiSlitModel from jwst.stpipe import query_step_status from jwst.assign_wcs import AssignWcsStep @@ -10,6 +10,9 @@ from jwst.extract_2d.tests.test_nirspec import create_nirspec_hdul from jwst.master_background import MasterBackgroundMosStep from jwst.master_background import nirspec_utils +from jwst.pixel_replace import PixelReplaceStep +from jwst.resample import ResampleSpecStep +from jwst.extract_1d import Extract1dStep def create_msa_hdul(): @@ -74,6 +77,7 @@ def nirspec_msa_metfl(tmp_path): hdul.close() return filename + @pytest.fixture def nirspec_msa_extracted2d(nirspec_msa_rate, nirspec_msa_metfl): model = ImageModel(nirspec_msa_rate) @@ -82,6 +86,20 @@ def nirspec_msa_extracted2d(nirspec_msa_rate, nirspec_msa_metfl): return model +def mk_multispec(model): + specs_model = MultiSlitModel() + specs_model.update(model) + slits = [] + for slit in model.slits: + if nirspec_utils.is_background_msa_slit(slit): + slits.append(slit) + specs_model.slits.extend(slits) + specs_model = PixelReplaceStep.call(specs_model) + specs_model = ResampleSpecStep.call(specs_model) + specs_model = Extract1dStep.call(specs_model) + return specs_model + + def test_master_background_mos(nirspec_msa_extracted2d): model = nirspec_msa_extracted2d @@ -97,43 +115,47 @@ def test_master_background_mos(nirspec_msa_extracted2d): # Check that a background was subtracted from the science data assert not np.allclose(sci_orig, sci_bkgsub) - model.close() - result.close() + del model + del result -def test_create_background_from_multislit(nirspec_msa_extracted2d): +def test_create_background_from_multispec(nirspec_msa_extracted2d): model = nirspec_msa_extracted2d - # Insert a outliers into one of the background spectra + # Insert outliers into one of the background spectra nypix = len(model.slits[0].data) nxpix = len(model.slits[0].data) - model.slits[0].data[nypix//2,nxpix//2-1:nxpix//2+1] = 10 + model.slits[0].data[nypix//2, nxpix//2-1:nxpix//2+1] = 10 + + specs_model = mk_multispec(model) # First check that we can make a master background from the inputs # Check that with sigma_clip=None, the outlier is retained - master_background, _ = nirspec_utils.create_background_from_multislit( - model, sigma_clip=None) + master_background = nirspec_utils.create_background_from_multispec( + specs_model, sigma_clip=None) assert np.any(master_background.spec[0].spec_table['surf_bright'] > 1) # Confirm that using a median_filter will filter out the outlier - master_background, _ = nirspec_utils.create_background_from_multislit( - model, median_kernel=4) + master_background = nirspec_utils.create_background_from_multispec( + specs_model, median_kernel=4) assert np.allclose(master_background.spec[0].spec_table['surf_bright'], 1) # Confirm that using a sigma clipping when combining background spectra # removes the outlier - master_background, _ = nirspec_utils.create_background_from_multislit( - model, sigma_clip=3) + master_background = nirspec_utils.create_background_from_multispec( + specs_model, sigma_clip=3) assert np.allclose(master_background.spec[0].spec_table['surf_bright'], 1) - model.close() + del model + del specs_model + def test_map_to_science_slits(nirspec_msa_extracted2d): model = nirspec_msa_extracted2d + specs_model = mk_multispec(model) - master_background, _ = nirspec_utils.create_background_from_multislit( - model) + master_background = nirspec_utils.create_background_from_multispec(specs_model) # Check that the master background is expanded to the shape of the input slits mb_multislit = nirspec_utils.map_to_science_slits(model, master_background) @@ -145,13 +167,15 @@ def test_map_to_science_slits(nirspec_msa_extracted2d): nonzero = slit_data != 0 assert np.allclose(slit_data[nonzero], 1) - model.close() + del model + del specs_model + def test_apply_master_background(nirspec_msa_extracted2d): model = nirspec_msa_extracted2d + specs_model = mk_multispec(model) - master_background, _ = nirspec_utils.create_background_from_multislit( - model) + master_background = nirspec_utils.create_background_from_multispec(specs_model) mb_multislit = nirspec_utils.map_to_science_slits(model, master_background) result = nirspec_utils.apply_master_background(model, mb_multislit, inverse=False) @@ -175,7 +199,8 @@ def test_apply_master_background(nirspec_msa_extracted2d): assert np.any(diff != 0) assert np.allclose(diff[diff != 0], -1) - model.close() - result.close() + del model + del result + del specs_model diff --git a/jwst/regtest/test_miri_mrs_spec3.py b/jwst/regtest/test_miri_mrs_spec3.py index 1b71a4c984..6e27960713 100644 --- a/jwst/regtest/test_miri_mrs_spec3.py +++ b/jwst/regtest/test_miri_mrs_spec3.py @@ -100,10 +100,14 @@ def run_spec3_ifushort_extract1d(rtdata_module): @pytest.mark.parametrize( 'output', [ - 'jw01024-c1000_t002_miri_ch3-mediumlong_x1d.fits', - 'jw01024-c1000_t002_miri_ch4-mediumlong_x1d.fits', - 'jw01024-c1000_t002_miri_ch3-mediumlong_s3d.fits', - 'jw01024-c1000_t002_miri_ch4-mediumlong_s3d.fits', + 'jw01024-c1000_t002_miri_ch3-long_x1d.fits', + 'jw01024-c1000_t002_miri_ch3-medium_x1d.fits', + 'jw01024-c1000_t002_miri_ch4-long_x1d.fits', + 'jw01024-c1000_t002_miri_ch4-medium_x1d.fits', + 'jw01024-c1000_t002_miri_ch3-long_s3d.fits', + 'jw01024-c1000_t002_miri_ch3-medium_s3d.fits', + 'jw01024-c1000_t002_miri_ch4-long_s3d.fits', + 'jw01024-c1000_t002_miri_ch4-medium_s3d.fits', ] ) def test_spec3_ifulong(run_spec3_ifulong, fitsdiff_default_kwargs, output): @@ -122,10 +126,14 @@ def test_spec3_ifulong(run_spec3_ifulong, fitsdiff_default_kwargs, output): @pytest.mark.parametrize( 'output', [ - 'jw01024-c1000_t002_miri_ch1-mediumlong_x1d.fits', - 'jw01024-c1000_t002_miri_ch2-mediumlong_x1d.fits', - 'jw01024-c1000_t002_miri_ch1-mediumlong_s3d.fits', - 'jw01024-c1000_t002_miri_ch2-mediumlong_s3d.fits' + 'jw01024-c1000_t002_miri_ch1-long_x1d.fits', + 'jw01024-c1000_t002_miri_ch1-medium_x1d.fits', + 'jw01024-c1000_t002_miri_ch2-long_x1d.fits', + 'jw01024-c1000_t002_miri_ch2-medium_x1d.fits', + 'jw01024-c1000_t002_miri_ch1-long_s3d.fits', + 'jw01024-c1000_t002_miri_ch1-medium_s3d.fits', + 'jw01024-c1000_t002_miri_ch2-long_s3d.fits', + 'jw01024-c1000_t002_miri_ch2-medium_s3d.fits' ], ) def test_spec3_ifushort(run_spec3_ifushort, fitsdiff_default_kwargs, output): @@ -145,10 +153,14 @@ def test_spec3_ifushort(run_spec3_ifushort, fitsdiff_default_kwargs, output): @pytest.mark.parametrize( 'output', [ - 'miri_mrs_emsm_ch1-mediumlong_x1d.fits', - 'miri_mrs_emsm_ch2-mediumlong_x1d.fits', - 'miri_mrs_emsm_ch1-mediumlong_s3d.fits', - 'miri_mrs_emsm_ch2-mediumlong_s3d.fits' + 'miri_mrs_emsm_ch1-long_x1d.fits', + 'miri_mrs_emsm_ch1-medium_x1d.fits', + 'miri_mrs_emsm_ch2-long_x1d.fits', + 'miri_mrs_emsm_ch2-medium_x1d.fits', + 'miri_mrs_emsm_ch1-long_s3d.fits', + 'miri_mrs_emsm_ch1-medium_s3d.fits', + 'miri_mrs_emsm_ch2-long_s3d.fits', + 'miri_mrs_emsm_ch2-medium_s3d.fits' ], ) def test_spec3_ifushort_emsm(run_spec3_ifushort_emsm, fitsdiff_default_kwargs, output): @@ -168,8 +180,10 @@ def test_spec3_ifushort_emsm(run_spec3_ifushort_emsm, fitsdiff_default_kwargs, o @pytest.mark.parametrize( 'output', [ - 'jw01024-c1000_t002_extract1dtest_miri_ch1-mediumlong_x1d.fits', - 'jw01024-c1000_t002_extract1dtest_miri_ch2-mediumlong_x1d.fits' + 'jw01024-c1000_t002_extract1dtest_miri_ch1-long_x1d.fits', + 'jw01024-c1000_t002_extract1dtest_miri_ch1-medium_x1d.fits', + 'jw01024-c1000_t002_extract1dtest_miri_ch2-long_x1d.fits', + 'jw01024-c1000_t002_extract1dtest_miri_ch2-medium_x1d.fits' ], ) def test_spec3_ifushort_extract1d(run_spec3_ifushort_extract1d, fitsdiff_default_kwargs, output): diff --git a/jwst/scripts/okify_regtests.py b/jwst/scripts/okify_regtests.py index 4a2494bb1d..4784ad1662 100755 --- a/jwst/scripts/okify_regtests.py +++ b/jwst/scripts/okify_regtests.py @@ -13,7 +13,6 @@ import tempfile from argparse import ArgumentParser from contextlib import contextmanager -from glob import glob from pathlib import Path import asdf @@ -110,11 +109,10 @@ def artifactory_get_breadcrumbs(build_number, suffix): args = list( ['jfrog', 'rt', 'dl'] + [f"{ARTIFACTORY_REPO}/*_GITHUB_CI_*-{build_number}/*{suffix}"] - + ['--flat'] ) subprocess.run(args, check=True, capture_output=True) - return sorted(glob(f'*{suffix}')) + return sorted([str(p) for p in Path().rglob(f'*{suffix}')]) def artifactory_get_build_artifacts(build_number): diff --git a/jwst/scripts/tests/test_scripts.py b/jwst/scripts/tests/test_scripts.py index b9ba4e8d08..eb133b15dc 100644 --- a/jwst/scripts/tests/test_scripts.py +++ b/jwst/scripts/tests/test_scripts.py @@ -9,7 +9,6 @@ 'asn_make_pool', 'collect_pipeline_cfgs', 'create_data', - 'okify_regtests', 'pointing_summary', 'schemadoc', 'set_telescope_pointing', diff --git a/pyproject.toml b/pyproject.toml index 26d333cffd..ef138c0999 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -18,7 +18,7 @@ classifiers = [ "Programming Language :: Python :: 3.12", ] dependencies = [ - "asdf>=3.1.0,<4", + "asdf>=3.1.0,<5", "astropy>=5.3", "BayesicFitting>=3.0.1", "crds>=12.0.3", @@ -74,7 +74,6 @@ collect_pipeline_cfgs = "jwst.scripts.collect_pipeline_cfgs:main" create_data = "jwst.scripts.create_data:main" csvconvert = "jwst.csv_tools.csvconvert:CSVConvertScript" exp_to_source = "jwst.exp_to_source.main:Main" -okify_regtests = "jwst.scripts.okify_regtests:main" pointing_summary = "jwst.scripts.pointing_summary:main" schemadoc = "jwst.scripts.schemadoc:main" set_telescope_pointing = "jwst.scripts.set_telescope_pointing:main" @@ -106,8 +105,6 @@ sdp = [ ] test = [ "ci-watson>=0.5.0", - "colorama>=0.4.1", - "readchar>=3.0", "ruff", "pysiaf>=0.13.0", "pytest>=6.0.0",