JP-3718: Improve AMI observable averaging, error calculations (#8846)

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.

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 <https://asdf-standard.readthedocs.io/>`_ 
+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 <https://doi.org/10.1051/0004-6361/201526405>`_.
 
-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".

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 <ami_analyze_step>` 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 <ami_normalize_step>` 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 <ami_analyze_step>` separately for access to flexible input arguments and to save additional diagnostic output products. See the step documentation for more details.
+

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}')

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  <rotate coords clockwise>
                 ctrs_rot = utils.rotate2dccw(mask_ctrs, np.deg2rad(-rot_ang))