JP-1737: Propagate error through Extract1D step

CHANGES.rst

@@ -106,6 +106,12 @@ documentation
 - Update ``assign_mtwcs`` step docs and reference the ``assign_mtwcs`` step in the
   ``calwebb_image3`` and ``calwebb_spec3`` pipeline docs [#6024]
 
+extract_1d
+----------
+
+- Implemented error and variance propagation for all modes but those
+  utilizing IFU cubes [#6014]
+
 extract_2d
 ----------
 

docs/jwst/data_products/science_products.rst

@@ -526,27 +526,45 @@ extension for each integration in the exposure.
 
 The structure of the "EXTRACT1D" table extension is as follows:
 
-+-------------+-----------+--------------------+---------------+
-| Column Name | Data Type | Contents           | Units         |
-+=============+===========+===================+================+
-| WAVELENGTH  | float64   | Wavelength values  | :math:`\mu` m |
-+-------------+-----------+--------------------+---------------+
-| FLUX        | float64   | Flux values        | Jy            |
-+-------------+-----------+--------------------+---------------+
-| ERROR       | float64   | Error values       | Jy            |
-+-------------+-----------+--------------------+---------------+
-| SURF_BRIGHT | float64   | Surface Brightness | MJy/sr        |
-+-------------+-----------+--------------------+---------------+
-| SB_ERROR    | float64   | Surf. Brt. errors  | MJy/sr        |
-+-------------+-----------+--------------------+---------------+
-| DQ          | uint32    | DQ flags           | N/A           |
-+-------------+-----------+--------------------+---------------+
-| BACKGROUND  | float64   | Background signal  | MJy/sr        |
-+-------------+-----------+--------------------+---------------+
-| BERROR      | float64   | Background error   | MJy/sr        |
-+-------------+-----------+--------------------+---------------+
-| NPIXELS     | float64   | Number of pixels   | N/A           |
-+-------------+-----------+--------------------+---------------+
++-------------------+-----------+--------------------+---------------+
+| Column Name       | Data Type | Contents           | Units         |
++===================+===========+===================+================+
+| WAVELENGTH        | float64   | Wavelength values  | :math:`\mu` m |
++-------------------+-----------+--------------------+---------------+
+| FLUX              | float64   | Flux values        | Jy            |
++-------------------+-----------+--------------------+---------------+
+| FLUX_ERROR        | float64   | Error values       | Jy            |
++-------------------+-----------+--------------------+---------------+
+| FLUX_VAR_POISSON  | float64   | Error values       | Jy^2          |
++-------------------+-----------+--------------------+---------------+
+| FLUX_VAR_RNOISE   | float64   | Error values       | Jy^2          |
++-------------------+-----------+--------------------+---------------+
+| FLUX_VAR_FLAT     | float64   | Error values       | Jy^2          |
++-------------------+-----------+--------------------+---------------+
+| SURF_BRIGHT       | float64   | Surface Brightness | MJy/sr        |
++-------------------+-----------+--------------------+---------------+
+| SB_ERROR          | float64   | Surf. Brt. errors  | MJy/sr        |
++-------------------+-----------+--------------------+---------------+
+| SB_VAR_POISSON    | float64   | Surf. Brt. errors  | (MJy/sr)^2    |
++-------------------+-----------+--------------------+---------------+
+| SB_VAR_RNOISE     | float64   | Surf. Brt. errors  | (MJy/sr)^2    |
++-------------------+-----------+--------------------+---------------+
+| SB_VAR_FLAT       | float64   | Surf. Brt. errors  | (MJy/sr)^2    |
++-------------------+-----------+--------------------+---------------+
+| DQ                | uint32    | DQ flags           | N/A           |
++-------------------+-----------+--------------------+---------------+
+| BACKGROUND        | float64   | Background signal  | MJy/sr        |
++-------------------+-----------+--------------------+---------------+
+| BKGD_ERROR        | float64   | Background error   | MJy/sr        |
++-------------------+-----------+--------------------+---------------+
+| BKGD_VAR_POISSON  | float64   | Background error   | (MJy/sr)^2    |
++-------------------+-----------+--------------------+---------------+
+| BKGD_VAR_RNOISE   | float64   | Background error   | (MJy/sr)^2    |
++-------------------+-----------+--------------------+---------------+
+| BKGD_VAR_FLAT     | float64   | Background error   | (MJy/sr)^2    |
++-------------------+-----------+--------------------+---------------+
+| NPIXELS           | float64   | Number of pixels   | N/A           |
++-------------------+-----------+--------------------+---------------+
 
 The table is constructed using a simple 2-D layout, using one row per extracted spectral
 element in the dispersion direction of the data (i.e. one row per wavelength bin).

docs/jwst/extract_1d/description.rst

@@ -37,8 +37,8 @@ results onto an infinite aperture scale.
 This is done by creating interpolation functions based on the APCORR reference
 file data and applying the interpolated aperture correction (a multiplicative
 factor between 0 and 1) to the extracted, 1D spectral data (corrected data
-include the "flux", "surf_bright", "error", and "sb_error" columns in the output
-table).
+include the "flux", "surf_bright", "flux_error", "sb_error", and all flux and
+surface brightness variance columns in the output table).
 
 Input
 -----
@@ -65,24 +65,27 @@ Output
 ------
 The output will be in ``MultiSpecModel`` format. For each input slit there will
 be an output table extension with the name EXTRACT1D.  This extension will
-have columns WAVELENGTH, FLUX, ERROR, SURF_BRIGHT, SB_ERROR, DQ,
-BACKGROUND, BERROR and NPIXELS.
+have columns WAVELENGTH, FLUX, FLUX_ERROR, FLUX_VAR_POISSON, FLUX_VAR_RNOISE,
+FLUX_VAR_FLAT, SURF_BRIGHT, SB_ERROR, SB_VAR_POISSON, SB_VAR_RNOISE,
+SB_VAR_FLAT, DQ, BACKGROUND, BKGD_ERROR, BKGD_VAR_POISSON, BKGD_VAR_RNOISE,
+BKGD_VAR_FLAT and NPIXELS.
 Some metadata will be written to the table header, mostly copied from the
 input header.
 
 The output WAVELENGTH data is copied from the wavelength array of the input 2D data,
 if that attribute exists and was populated, otherwise it is calculated from the WCS.
 FLUX is the flux density in Janskys; see keyword TUNIT2 if the data are
-in a FITS BINTABLE.  ERROR is the error estimate for FLUX, and it has the
-same units as FLUX.
+in a FITS BINTABLE.  FLUX_ERROR is the error estimate for FLUX, and it has the
+same units as FLUX. The error is calculated as the square root of the sum of the
+three variance arrays: Poisson, read noise (RNOISE), and flat field (FLAT).
 SURF_BRIGHT is the surface brightness in MJy / sr, except that for point
 sources observed with NIRSpec and NIRISS SOSS, SURF_BRIGHT will be set to
 zero, because there's no way to express the extracted results from those modes
-as a surface brightness. SB_ERROR is the error estimate for SURF_BRIGHT.
-While it's expected that a user will make use of the FLUX column for
-point-source data and the SURF_BRIGHT column for an extended source,
-both columns are populated (except for NIRSpec and NIRISS SOSS point sources,
-as mentioned above).
+as a surface brightness. SB_ERROR is the error estimate for SURF_BRIGHT, calculated
+in the same fashion as FLUX_ERROR but using the SB_VAR arrays. While it's expected
+that a user will make use of the FLUX column for point-source data and the
+SURF_BRIGHT column for an extended source, both columns are populated (except for
+NIRSpec and NIRISS SOSS point sources, as mentioned above).
 The ``extract_1d`` step collapses the input data from 2-D to 1-D by summing
 one or more rows (or columns, depending on the dispersion direction).
 A background may optionally be subtracted, but
@@ -104,8 +107,8 @@ extraction region.  Note that this is not necessarily a constant, and
 the value is not necessarily an integer (the data type is float).
 BACKGROUND is the measured background, scaled to the extraction width used
 for FLUX and SURF_BRIGHT.  BACKGROUND will be zero if background subtraction
-is not requested.
-ERROR, SB_ERROR, BERROR, and DQ are not populated with useful values yet.
+is not requested. BKGD_ERROR is calculated as the square root of the sum of the
+BKGD_VAR arrays. DQ is not populated with useful values yet.
 
 
 .. _extract-1d-for-slits:

jwst/combine_1d/combine1d.py

@@ -17,7 +17,7 @@ class InputSpectrumModel:
     """Attributes:
         wavelength
         flux
-        error
+        flux_error
         surf_bright
         sb_error
         dq
@@ -51,7 +51,7 @@ def __init__(self, ms, spec, exptime_key):
         self.wavelength = spec.spec_table.field("wavelength").copy()
 
         self.flux = spec.spec_table.field("flux").copy()
-        self.error = spec.spec_table.field("error").copy()
+        self.flux_error = spec.spec_table.field("flux_error").copy()
         try:
             self.surf_bright = spec.spec_table.field("surf_bright").copy()
             self.sb_error = spec.spec_table.field("sb_error").copy()
@@ -87,7 +87,7 @@ def __init__(self, ms, spec, exptime_key):
     def close(self):
         self.wavelength = None
         self.flux = None
-        self.error = None
+        self.flux_error = None
         self.surf_bright = None
         self.sb_error = None
         self.dq = None
@@ -103,7 +103,7 @@ class OutputSpectrumModel:
     """Attributes:
         wavelength
         flux
-        error
+        flux_error
         surf_bright
         sb_error
         dq
@@ -117,7 +117,7 @@ def __init__(self):
 
         self.wavelength = None
         self.flux = None
-        self.error = None
+        self.flux_error = None
         self.surf_bright = None
         self.sb_error = None
         self.dq = None
@@ -181,7 +181,7 @@ def accumulate_sums(self, input_spectra):
         nelem = self.wavelength.shape[0]
 
         self.flux = np.zeros(nelem, dtype=np.float64)
-        self.error = np.zeros(nelem, dtype=np.float64)
+        self.flux_error = np.zeros(nelem, dtype=np.float64)
         self.surf_bright = np.zeros(nelem, dtype=np.float64)
         self.sb_error = np.zeros(nelem, dtype=np.float64)
         self.dq = np.zeros(nelem, dtype=dq_dtype)
@@ -214,7 +214,7 @@ def accumulate_sums(self, input_spectra):
                 weight = in_spec.weight[i]
                 self.dq[k] |= in_spec.dq[i]
                 self.flux[k] += in_spec.flux[i] * weight
-                self.error[k] += (in_spec.error[i] * weight)**2
+                self.flux_error[k] += (in_spec.flux_error[i] * weight)**2
                 self.surf_bright[k] += (in_spec.surf_bright[i] * weight)
                 self.sb_error[k] += (in_spec.sb_error[i] * weight)**2
                 self.weight[k] += weight
@@ -234,7 +234,7 @@ def accumulate_sums(self, input_spectra):
             log.warning("    these elements will be omitted.")
             self.wavelength = self.wavelength[index]
             self.flux = self.flux[index]
-            self.error = self.error[index]
+            self.flux_error = self.flux_error[index]
             self.surf_bright = self.surf_bright[index]
             self.sb_error = self.sb_error[index]
             self.dq = self.dq[index]
@@ -251,7 +251,7 @@ def compute_combination(self):
             sum_weight = np.where(self.weight > 0., self.weight, 1.)
             self.surf_bright /= sum_weight
             self.flux /= sum_weight
-            self.error = np.sqrt(self.error / sum_weight)
+            self.flux_error = np.sqrt(self.flux_error / sum_weight)
             self.sb_error = np.sqrt(self.sb_error / sum_weight)
             self.normalized = True
 
@@ -273,7 +273,7 @@ def create_output_data(self):
         # Note that these arrays have to be in the right order.
         data = np.array(list(zip(self.wavelength,
                                  self.flux,
-                                 self.error,
+                                 self.flux_error,
                                  self.surf_bright,
                                  self.sb_error,
                                  self.dq,
@@ -286,7 +286,7 @@ def create_output_data(self):
     def close(self):
         self.wavelength = None
         self.flux = None
-        self.error = None
+        self.flux_error = None
         self.surf_bright = None
         self.sb_error = None
         self.dq = None

jwst/combine_1d/tests/test_dq.py

@@ -38,20 +38,25 @@ def create_spec_model(npoints=10, flux=1e-9, wave_range=(11, 13)):
     error = np.zeros(npoints)
     surf_bright = np.zeros(npoints)
     sb_error = np.zeros(npoints)
+    var_dummy = error.copy()
     dq = np.zeros(npoints)
     background = np.zeros(npoints)
     berror = np.zeros(npoints)
     npixels = np.zeros(npoints)
 
-    data = [(i, j, k, l, m, n, o, p, q) for i, j, k, l, m, n, o, p, q in zip(wavelength, flux, error,
-                                                                             surf_bright, sb_error, dq,
-                                                                             background, berror, npixels)]
+    spec_dtype = datamodels.SpecModel().spec_table.dtype  # This data type is used for creating an output table.
 
-    spec_table = np.array(data, dtype=[('WAVELENGTH', 'f8'), ('FLUX', 'f8'),
-                                       ('ERROR', 'f8'), ('SURF_BRIGHT', 'f8'),
-                                       ('SB_ERROR', 'f8'), ('DQ', 'u4'), ('BACKGROUND', 'f8'),
-                                       ('BERROR', 'f8'), ('NPIXELS', 'f8')])
-    spec_model = datamodels.SpecModel()
-    spec_model.spec_table = spec_table
+    otab = np.array(
+        list(
+            zip(
+                wavelength, flux, error, var_dummy, var_dummy, var_dummy,
+                surf_bright, sb_error, var_dummy, var_dummy, var_dummy,
+                dq, background, berror, var_dummy, var_dummy, var_dummy,
+                npixels
+            ),
+        ), dtype=spec_dtype
+    )
+
+    spec_model = datamodels.SpecModel(spec_table=otab)
 
     return spec_model

jwst/datamodels/schemas/spectable.schema.yaml

@@ -9,17 +9,35 @@ datatype:
   datatype: float64
 - name: FLUX
   datatype: float64
-- name: ERROR
+- name: FLUX_ERROR
+  datatype: float64
+- name: FLUX_VAR_POISSON
+  datatype: float64
+- name: FLUX_VAR_RNOISE
+  datatype: float64
+- name: FLUX_VAR_FLAT
   datatype: float64
 - name: SURF_BRIGHT
   datatype: float64
 - name: SB_ERROR
   datatype: float64
+- name: SB_VAR_POISSON
+  datatype: float64
+- name: SB_VAR_RNOISE
+  datatype: float64
+- name: SB_VAR_FLAT
+  datatype: float64
 - name: DQ
   datatype: uint32
 - name: BACKGROUND
   datatype: float64
-- name: BERROR
+- name: BKGD_ERROR
+  datatype: float64
+- name: BKGD_VAR_POISSON
+  datatype: float64
+- name: BKGD_VAR_RNOISE
+  datatype: float64
+- name: BKGD_VAR_FLAT
   datatype: float64
 - name: NPIXELS
   datatype: float64

jwst/extract_1d/apply_apcorr.py

@@ -153,7 +153,9 @@ def apply(self, spec_table: fits.FITS_rec):
             Table of aperture corrections values from apcorr reference file.
 
         """
-        cols_to_correct = ('flux', 'surf_bright', 'error', 'sb_error')
+        cols_to_correct = ('flux', 'flux_error', 'flux_var_poisson', 'flux_var_rnoise',
+                           'flux_var_flat', 'surf_bright', 'sb_error', 'sb_var_poisson',
+                           'sb_var_rnoise', 'sb_var_flat')
 
         for row in spec_table:
             correction = self.apcorr_func(row['npixels'], row['wavelength'])
@@ -227,7 +229,9 @@ def apply(self, spec_table: fits.FITS_rec):
             Table of aperture corrections values from apcorr reference file.
 
         """
-        cols_to_correct = ('flux', 'surf_bright', 'error', 'sb_error')
+        cols_to_correct = ('flux', 'flux_error', 'flux_var_poisson', 'flux_var_rnoise',
+                           'flux_var_flat', 'surf_bright', 'sb_error', 'sb_var_poisson',
+                           'sb_var_rnoise', 'sb_var_flat')
 
         for row in spec_table:
             try:
@@ -284,7 +288,9 @@ def apply(self, spec_table: fits.FITS_rec):
             Table of aperture corrections values from apcorr reference file.
 
         """
-        cols_to_correct = ('flux', 'surf_bright', 'error', 'sb_error')
+        cols_to_correct = ('flux', 'flux_error', 'flux_var_poisson', 'flux_var_rnoise',
+                           'flux_var_flat', 'surf_bright', 'sb_error', 'sb_var_poisson',
+                           'sb_var_rnoise', 'sb_var_flat')
 
         for i, row in enumerate(spec_table):
             correction = self.apcorr_correction[i]