JP-3901: Enable Detector1Pipeline processing of DHS data

jwst/jump/jump_step.py

@@ -158,13 +158,13 @@
                 gain_2d = gain_m.data
             else:
                 self.log.info("Extracting gain subarray to match science data")
-                gain_2d = reffile_utils.get_subarray_data(result, gain_m)
+                gain_2d = reffile_utils.get_subarray_model(result, gain_m).data
 
             if reffile_utils.ref_matches_sci(result, rnoise_m):
                 rnoise_2d = rnoise_m.data
             else:
                 self.log.info("Extracting readnoise subarray to match science data")
-                rnoise_2d = reffile_utils.get_subarray_data(result, rnoise_m)
+                rnoise_2d = reffile_utils.get_subarray_model(result, rnoise_m).data
 
         # Instantiate a JumpData class and populate it based on the input RampModel.
         jump_data = JumpData(result, gain_2d, rnoise_2d, dqflags.pixel)

jwst/lib/reffile_utils.py

@@ -197,101 +197,6 @@
         return False
 
 
-def get_subarray_data(sci_model, ref_model):
-    """
-    Extract a subarray from the data attribute of a reference file
-    data model that matches the subarray characteristics of a
-    science data model. Only the extracted data array is returned.
-
-    Parameters
-    ----------
-    sci_model: JWST data model
-        science data model
-
-    ref_model: JWST data model
-        reference file data model
-
-    Returns
-    -------
-    array: 2-D extracted data array
-    """
-
-    # Make sure xstart/ystart exist in science data model
-    if sci_model.meta.subarray.xstart is None or sci_model.meta.subarray.ystart is None:
-
-        # If the science file is full-frame, set the missing params to
-        # default values
-        if sci_model.data.shape[-1] == 2048 and sci_model.data.shape[-2] == 2048:
-            sci_model.meta.subarray.xstart = 1
-            sci_model.meta.subarray.ystart = 1
-            sci_model.meta.subarray.xsize = 2048
-            sci_model.meta.subarray.ysize = 2048
-        else:
-            raise ValueError('xstart or ystart metadata values not found in input model')
-
-    # Make sure xstart/ystart exist in reference data model
-    if ref_model.meta.subarray.xstart is None or ref_model.meta.subarray.ystart is None:
-
-        # If the ref file is full-frame, set the missing params to
-        # default values
-        if ref_model.meta.instrument.name.upper() == 'MIRI':
-            if ref_model.data.shape[-1] == 1032 and ref_model.data.shape[-2] == 1024:
-                ref_model.meta.subarray.xstart = 1
-                ref_model.meta.subarray.ystart = 1
-                ref_model.meta.subarray.xsize = 1032
-                ref_model.meta.subarray.ysize = 1024
-            else:
-                raise ValueError('xstart or ystart metadata values not found in reference model')
-        else:
-            if ref_model.data.shape[-1] == 2048 and ref_model.data.shape[-2] == 2048:
-                ref_model.meta.subarray.xstart = 1
-                ref_model.meta.subarray.ystart = 1
-                ref_model.meta.subarray.xsize = 2048
-                ref_model.meta.subarray.ysize = 2048
-            else:
-                raise ValueError('xstart or ystart metadata values not found in reference model')
-
-    # Get subarray limits from metadata of input model
-    xstart_sci = sci_model.meta.subarray.xstart
-    xsize_sci = sci_model.meta.subarray.xsize
-    ystart_sci = sci_model.meta.subarray.ystart
-    ysize_sci = sci_model.meta.subarray.ysize
-    log.debug('science xstart=%d, xsize=%d, ystart=%d, ysize=%d',
-              xstart_sci, xsize_sci, ystart_sci, ysize_sci)
-
-    # Get subarray limits from metadata of reference model
-    xstart_ref = ref_model.meta.subarray.xstart
-    xsize_ref = ref_model.meta.subarray.xsize
-    ystart_ref = ref_model.meta.subarray.ystart
-    ysize_ref = ref_model.meta.subarray.ysize
-    log.debug('reference xstart=%d, xsize=%d, ystart=%d, ysize=%d',
-              xstart_ref, xsize_ref, ystart_ref, ysize_ref)
-
-    # Compute slice limits, in 0-indexed python notation
-    xstart = xstart_sci - xstart_ref
-    ystart = ystart_sci - ystart_ref
-    xstop = xstart + xsize_sci
-    ystop = ystart + ysize_sci
-    log.debug('slice xstart=%d, xstop=%d, ystart=%d, ystop=%d',
-              xstart, xstop, ystart, ystop)
-
-    # Make sure that the slice limits are within the bounds of
-    # the reference file data array
-    if (xstart < 0 or ystart < 0 or xstop > ref_model.meta.subarray.xsize or
-            ystop > ref_model.meta.subarray.ysize):
-        log.error('Computed reference file slice indexes are ' +
-                  'incompatible with size of reference data array')
-        log.error('Science: SUBSTRT1=%d, SUBSTRT2=%d, SUBSIZE1=%d, SUBSIZE2=%d',
-                  xstart_sci, ystart_sci, xsize_sci, ysize_sci)
-        log.error('Reference: SUBSTRT1=%d, SUBSTRT2=%d, SUBSIZE1=%d, SUBSIZE2=%d',
-                  xstart_ref, ystart_ref, xsize_ref, ysize_ref)
-        log.error('Slice indexes: xstart=%d, xstop=%d, ystart=%d, ystop=%d',
-                  xstart, xstop, ystart, ystop)
-        raise ValueError('Bad reference file slice indexes')
-
-    return ref_model.data[ystart:ystop, xstart:xstop]
-
-
 def get_subarray_model(sci_model, ref_model):
     """
     Create a subarray version of a reference file model that matches
@@ -312,6 +217,10 @@
     sub_model: JWST data model
         subarray version of the reference file model
     """
+    # If science data is in multistripe readout, use
+    # multistripe-specific subarray reconstruction.
+    if getattr(sci_model.meta.subarray, 'multistripe_reads1') is not None:
+        return get_multistripe_subarray_model(sci_model, ref_model)
 
     # Get the science model subarray params
     xstart_sci = sci_model.meta.subarray.xstart
@@ -388,6 +297,236 @@
     return sub_model
 
 
+def get_multistripe_subarray_model(sci_model, ref_model):
+    """
+    Create a subarray version of a reference file model that matches
+    the multistripe-specific subarray characteristics of a science
+    data model. A new model is created that contains subarrays of
+    all data arrays contained in the reference file model.
+
+    Parameters
+    ----------
+    sci_model: JWST data model
+        science data model
+
+    ref_model: JWST data model
+        reference file data model
+
+    Returns
+    -------
+    sub_model: JWST data model
+        subarray version of the reference file model
+    """
+    # Could consider just brute-forcing every attribute in a list,
+    # catching any errors for those attribs that don't exist for
+    # a given reftype? Catching unknown reftypes is probably safer,
+    # given possibility of unknown attrib.
+    if isinstance(ref_model, datamodels.MaskModel):
+        sub_model = stripe_read(sci_model, ref_model, ['dq'])
+    elif isinstance(ref_model, datamodels.GainModel):
+        sub_model = stripe_read(sci_model, ref_model, ['data'])
+    elif isinstance(ref_model, datamodels.LinearityModel):
+        sub_model = stripe_read(sci_model, ref_model, ['coeffs', 'dq'])
+    elif isinstance(ref_model, datamodels.ReadnoiseModel):
+        sub_model = stripe_read(sci_model, ref_model, ['data'])
+    elif isinstance(ref_model, datamodels.SaturationModel):
+        sub_model = stripe_read(sci_model, ref_model, ['data', 'dq'])
+    elif isinstance(ref_model, datamodels.SuperBiasModel):
+        sub_model = stripe_read(sci_model, ref_model, ['data', 'err', 'dq'])
+    else:
+        log.warning('Unsupported reference file model type')
+        sub_model = None
+
+    # Temporary saving of all sub_models created
+    sub_model.write(
+        sub_model.meta.filename.split('.')[0]
+        + '_'
+        + sci_model.meta.subarray.name
+        + '.fits'
+    )
+    return sub_model
+
+
+def stripe_read(sci_model, ref_model, attribs):
+    # Get the science model multistripe params
+    nreads1 = sci_model.meta.subarray.multistripe_reads1
+    nskips1 = sci_model.meta.subarray.multistripe_skips1
+    nreads2 = sci_model.meta.subarray.multistripe_reads2
+    nskips2 = sci_model.meta.subarray.multistripe_skips2
+    repeat_stripe = sci_model.meta.subarray.repeat_stripe
+    interleave_reads1 = sci_model.meta.subarray.interleave_reads1
+    xsize_sci = sci_model.meta.subarray.xsize
+    ysize_sci = sci_model.meta.subarray.ysize
+
+    # Get the reference model subarray params
+    sub_model = type(ref_model)()
+    sub_model.update(ref_model)
+    for attrib in attribs:
+        ref_array = getattr(ref_model, attrib)
+
+        sub_model[attrib] = generate_stripe_array(
+            ref_array,
+            xsize_sci,
+            ysize_sci,
+            nreads1,
+            nreads2,
+            nskips1,
+            nskips2,
+            repeat_stripe,
+            interleave_reads1,
+            sci_model.meta.subarray.fastaxis,
+            sci_model.meta.subarray.slowaxis,
+        )
+    return sub_model
+
+
+def generate_stripe_array(
+        ref_array,
+        xsize_sci,
+        ysize_sci,
+        nreads1,
+        nreads2,
+        nskips1,
+        nskips2,
+        repeat_stripe,
+        interleave_reads1,
+        fastaxis,
+        slowaxis,
+):
+    """
+
+    Parameters
+    ----------
+    ref_array : np.array
+        The scene to be sliced.
+    xsize_sci : int
+        Output shape in x dim.
+    ysize_sci : int
+        Output shape in y dim.
+    nreads1 : int
+        Multistripe header keyword.
+    nreads2 : int
+        Multistripe header keyword.
+    nskips1 : int
+        Multistripe header keyword.
+    nskips2 : int
+        Multistripe header keyword.
+    repeat_stripe : int
+        Multistripe header keyword.
+    interleave_reads1 : int
+        Multistripe header keyword.
+    fastaxis : int
+        The subarray keyword describing
+        the fast readout axis and direction.
+    slowaxis : int
+        The subarray keyword describing
+        the slow readout axis and direction.
+
+    Returns
+    -------
+
+    """
+    # Transform science data to detector frame
+    ref_array = science_detector_frame_transform(ref_array, fastaxis, slowaxis)
+    ref_shape = np.shape(ref_array)
+    stripe_out = np.zeros(
+        (*ref_shape[:-2], ysize_sci, xsize_sci), dtype=getattr(ref_array, 'dtype')
+    )
+    # Track the read position in the full frame with linecount, and number of lines
+    # read into subarray with sub_lines
+    linecount = 0
+    sub_lines = 0
+
+    # Start at 0, make nreads1 row reads
+    stripe_out[..., sub_lines: sub_lines + nreads1, :] = \
+        ref_array[..., linecount: linecount + nreads1, :]
+    linecount += nreads1
+    sub_lines += nreads1
+    # Now skip nskips1
+    linecount += nskips1
+    # Nreads2
+    stripe_out[..., sub_lines: sub_lines + nreads2, :] = \
+        ref_array[..., linecount: linecount + nreads2, :]
+    linecount += nreads2
+    sub_lines += nreads2
+
+    # Now, while the output size is less than the science array size:
+    # 1a. If repeat_stripe, reset linecount (HEAD) to initial position
+    #     after every nreads2.
+    # 1b. Else, do nskips2 followed by nreads2 until subarray complete.
+    # 2.  Following 1a., repeat sequence of nreads1, skips*, nreads2
+    #     until complete. For skips*:
+    # 3a. If interleave_reads1, value of skips increments by nskips2
+    #     for each stripe read.
+    # 3b. If not interleave, each loop after linecount reset is simply
+    #     nreads1 + nskips1 + nreads2.
+    interleave_skips = nskips1
+    while sub_lines < ysize_sci:
+        # If repeat_stripe, add interleaved rows to output and increment sub_lines
+        if repeat_stripe > 0:
+            linecount = 0
+            stripe_out[..., sub_lines: sub_lines + nreads1, :] = \
+                ref_array[..., linecount: linecount + nreads1, :]
+            linecount += nreads1
+            sub_lines += nreads1
+            if interleave_reads1:
+                interleave_skips += nskips2
+                linecount += interleave_skips
+            else:
+                linecount += nskips1
+        else:
+            linecount += nskips2
+        stripe_out[..., sub_lines: sub_lines + nreads2, :] = \
+            ref_array[..., linecount: linecount + nreads2, :]
+        linecount += nreads2
+        sub_lines += nreads2
+
+    if sub_lines != ysize_sci:
+        raise ValueError('Stripe readout resulted in mismatched reference array shape '
+                         'with respect to science array!')
+
+    # Transform from detector frame back to science frame
+    stripe_out = science_detector_frame_transform(stripe_out, fastaxis, slowaxis)
+
+    return stripe_out
+
+
+def science_detector_frame_transform(data, fastaxis, slowaxis):
+    """
+    Swap data array between science and detector frames.
+
+    Use the fastaxis and slowaxis keywords to invert
+    and/or transpose data array axes to move between the
+    science frame and the detector frame.
+
+    Parameters
+    ----------
+    data : np.array
+        Science array containing at least two dimensions.
+    fastaxis : int
+        Value of the fastaxis keyword for the data array
+        to be transformed.
+    slowaxis : int
+        Value of the slowaxis keyword for the data array
+        to be transformed.
+
+    Returns
+    -------
+    np.array
+        Data array transformed between science and
+        detector frames.
+    """
+    # If fastaxis is x-axis
+    if np.abs(fastaxis) == 1:
+        # Use sign of keywords to possibly reverse the ordering of the axes.
+        data = data[..., ::slowaxis//np.abs(slowaxis), ::fastaxis//np.abs(fastaxis)]
+    # Else fastaxis is y-axis, also need to transpose array
+    else:
+        data = data[..., ::fastaxis//np.abs(fastaxis), ::slowaxis//np.abs(slowaxis)]
+        data = np.swapaxes(data, -2, -1)
+    return data
+
+
 class MatchRowError(Exception):
     """
     Raised when more than one row is matched in a FITS table or list of dict.

jwst/ramp_fitting/ramp_fit_step.py

@@ -52,13 +52,13 @@ def get_reference_file_subarrays(model, readnoise_model, gain_model):
         gain_2d = gain_model.data
     else:
         log.info("Extracting gain subarray to match science data")
-        gain_2d = reffile_utils.get_subarray_data(model, gain_model)
+        gain_2d = reffile_utils.get_subarray_model(model, gain_model).data
 
     if reffile_utils.ref_matches_sci(model, readnoise_model):
         readnoise_2d = readnoise_model.data
     else:
         log.info("Extracting readnoise subarray to match science data")
-        readnoise_2d = reffile_utils.get_subarray_data(model, readnoise_model)
+        readnoise_2d = reffile_utils.get_subarray_model(model, readnoise_model).data
 
     return readnoise_2d, gain_2d