Merge branch 'main' into hackday_mteodoro_sky_variance_readnoise

changes/1585.mosaic_pipeline.rst

@@ -0,0 +1 @@
+Roundtrip L3 wcsinfo especially when skycell specifications are used

changes/1593.exposure_pipeline.rst

@@ -0,0 +1 @@
+Fix description of arguments in docs and add description of fully saturated input processing.

changes/1593.resample.rst

@@ -0,0 +1 @@
+Remove unused arguments from step specification.

changes/1593.saturation.rst

@@ -0,0 +1 @@
+Add saturation step docs to package index.

docs/roman/package_index.rst

@@ -14,6 +14,7 @@
    references_general/index.rst
    refpix/index.rst
    resample/index.rst
+   saturation/index.rst
    skymatch/index.rst
    source_catalog/index.rst
    stpipe/index.rst

docs/roman/pipeline/exposure_pipeline.rst

@@ -39,18 +39,7 @@ table below.
 
 Arguments
 ---------
-The ``exposure`` pipeline has an optional argument::
-
-  --use_ramp_jump_detection  boolean  default=True
-
-When set to ``True``, the pipeline will perform jump detection as a part of the ramp
-fitting  step. The data at this stage of the pipeline are still in the form of the original
-3D ramps ( ngroups x ncols x nrows ) and have had all of the detector-level
-correction steps applied to it, up to but not including the detection and flagging of
-Cosmic-Ray (CR) hits within each ramp (integration). For this case the  jump detection
-module in :ref:`ramp_fitting <ramp_fitting_step>` will update the dq array with the CR hits (jumps) that
-are identified in the step.
-
+The ``exposure`` pipeline has no arguments
 
 Inputs
 ------
@@ -78,6 +67,16 @@ extensions. When such a file is loaded into the pipeline, it is immediately
 converted into a `~romancal.datamodels.RampModel`, and has all additional data arrays
 for errors and Data Quality flags created and initialized.
 
+When the ``ExposurePipeline`` processes a fully saturated input (all pixels flagged as saturated).
+The corresponding output image will:
+
+- contain all 0 data arrays
+- contain all 0 variance arrays
+- not be processed by steps beyond saturation
+
+A single fully saturated input will also cause :ref:`tweakreg <tweakreg_step>` to be skipped
+for all input images.
+
 Outputs
 -------
 

docs/roman/references_general/references_general.rst

@@ -51,8 +51,6 @@ documentation on each reference file.
 +---------------------------------------------+--------------------------------------------------+
 | :ref:`flatfield <flatfield_step>`           | :ref:`FLAT <flat_reffile>`                       |
 +---------------------------------------------+--------------------------------------------------+
-| :ref:`GAIN <gain_reffile>`                       | :ref:`ramp_fitting <ramp_fitting_step>`     |
-+---------------------------------------------+--------------------------------------------------+
 | :ref:`linearity <linearity_step>`           | :ref:`LINEARITY <linearity_reffile>`             |
 +---------------------------------------------+--------------------------------------------------+
 | :ref:`photom <photom_step>`                 | :ref:`PHOTOM <photom_reffile>`                   |
@@ -74,6 +72,8 @@ documentation on each reference file.
 +--------------------------------------------------+---------------------------------------------+
 | :ref:`FLAT <flat_reffile>`                       | :ref:`flatfield <flatfield_step>`           |
 +--------------------------------------------------+---------------------------------------------+
+| :ref:`GAIN <gain_reffile>`                       | :ref:`ramp_fitting <ramp_fitting_step>`     |
++--------------------------------------------------+---------------------------------------------+
 | :ref:`LINEARITY <linearity_reffile>`             | :ref:`linearity <linearity_step>`           |
 +--------------------------------------------------+---------------------------------------------+
 | :ref:`MASK <mask_reffile>`                       | :ref:`dq_init <dq_init_step>`               |

docs/roman/resample/arguments.rst

@@ -92,24 +92,6 @@ image.
     where :math:`\text{DATA}` and :math:`med(\text{DATA})` correspond to the data array and its median, respectively.
     If the ``VAR_SKY`` array does not exist, the weight is set to 1 for all pixels (equal weighting).
 
-``--single`` (bool, default=False)
-    If set to `True`, resample each input image into a separate output.  If
-    `False` (the default), each input is resampled additively (with weights) to
-    a common output
-
-``--blendheaders`` (bool, default=True)
-    Blend metadata from all input images into the resampled output image.
-
-``--allowed_memory`` (float, default=None)
-    Specifies the fractional amount of free memory to allow when creating the
-    resampled image. If ``None``, the environment variable
-    ``DMODEL_ALLOWED_MEMORY`` is used. If not defined, no check is made. If the
-    resampled image would be larger than specified, an ``OutputTooLargeError``
-    exception will be generated.
-
-    For example, if set to ``0.5``, only resampled images that use less than
-    half the available memory can be created.
-
 ``--in_memory`` (bool, default=True)
     If set to `False`, write output datamodel to disk.
 

docs/roman/stpipe/config_asdf.rst

@@ -70,7 +70,6 @@ entries for clarity:
       output_use_model: false
       post_hooks: []
       pre_hooks: []
-      save_calibrated_ramp: false
       save_results: true
       search_output_file: true
       skip: false

docs/roman/stpipe/user_pipeline.rst

@@ -63,7 +63,6 @@ class:
     output_use_model: false
     post_hooks: []
     pre_hooks: []
-    save_calibrated_ramp: false
     save_results: true
     search_output_file: true
     skip: false

romancal/pipeline/exposure_pipeline.py

@@ -46,7 +46,6 @@ class ExposurePipeline(RomanPipeline):
     class_alias = "roman_elp"
 
     spec = """
-        save_calibrated_ramp = boolean(default=False)
         save_results = boolean(default=False)
         suffix = string(default="cal")
     """

romancal/pipeline/mosaic_pipeline.py

@@ -13,6 +13,7 @@
 from astropy import units as u
 from astropy.modeling import models
 from gwcs import WCS, coordinate_frames
+from stcal.alignment import util as wcs_util
 
 import romancal.datamodels.filetype as filetype
 from romancal.datamodels import ModelLibrary
@@ -114,20 +115,21 @@ def process(self, input):
 
                     # check to see if there exists a skycell on disk if not create it
                     if not isfile(skycell_file_name):
-                        # extract the wcs info from the record for generate_tan_wcs
+                        # extract the wcs info from the record for skycell_to_wcs
                         log.info(
                             "Creating skycell image at ra: %f  dec %f",
                             float(skycell_record["ra_center"]),
                             float(skycell_record["dec_center"]),
                         )
-                        skycell_wcs = generate_tan_wcs(skycell_record)
+                        skycell_wcs = skycell_to_wcs(skycell_record)
                         # skycell_wcs.bounding_box = bounding_box
 
                         # For resample to use an external grid we need to pass it the skycell gwcs object
                         # Currently we cannot do that directly so create an asdf file to read the skycell gwcs object
                         wcs_tree = {"wcs": skycell_wcs}
                         wcs_file = asdf.AsdfFile(wcs_tree)
                         wcs_file.write_to("skycell_wcs.asdf")
+
                         self.resample.output_wcs = "skycell_wcs.asdf"
                         self.resample.output_shape = (
                             int(skycell_record["nx"]),
@@ -162,41 +164,96 @@ def process(self, input):
         return result
 
 
-def generate_tan_wcs(skycell_record):
-    """extract the wcs info from the record for generate_tan_wcs
-    we need the scale, ra, dec, bounding_box"""
+def skycell_to_wcs(skycell_record):
+    """From a skycell record, generate a GWCS
+
+    Parameters
+    ----------
+    skycell_record : dict
+        A skycell record, or row, from the skycell patches table.
+
+    Returns
+    -------
+    wcsobj : wcs.GWCS
+        The GWCS object from the skycell record.
+    """
+    wcsinfo = dict()
+
+    # The scale is given in arcseconds per pixel. Convert to degrees.
+    wcsinfo["pixel_scale"] = float(skycell_record["pixel_scale"]) / 3600.0
+
+    # Remaining components of the wcsinfo block
+    wcsinfo["ra_ref"] = float(skycell_record["ra_projection_center"])
+    wcsinfo["dec_ref"] = float(skycell_record["dec_projection_center"])
+    wcsinfo["x_ref"] = float(skycell_record["x0_projection"])
+    wcsinfo["y_ref"] = float(skycell_record["y0_projection"])
+    wcsinfo["orientat"] = float(skycell_record["orientat_projection_center"])
+    wcsinfo["rotation_matrix"] = None
 
-    scale = float(skycell_record["pixel_scale"])
-    ra_center = float(skycell_record["ra_projection_center"])
-    dec_center = float(skycell_record["dec_projection_center"])
-    shiftx = float(skycell_record["x0_projection"])
-    shifty = float(skycell_record["y0_projection"])
+    # Bounding box of the skycell. Note that the center of the pixels are at (0.5, 0.5)
     bounding_box = (
         (-0.5, -0.5 + skycell_record["nx"]),
         (-0.5, -0.5 + skycell_record["ny"]),
     )
 
-    # components of the model
-    # shift = models.Shift(shiftx) & models.Shift(shifty)
+    wcsobj = wcsinfo_to_wcs(wcsinfo, bounding_box=bounding_box)
+    return wcsobj
+
+
+def wcsinfo_to_wcs(wcsinfo, bounding_box=None, name="wcsinfo"):
+    """Create a GWCS from the L3 wcsinfo meta
+
+    Parameters
+    ----------
+    wcsinfo : dict or MosaicModel.meta.wcsinfo
+        The L3 wcsinfo to create a GWCS from.
 
-    # select a scale for the skycell image, this will come from INS and may
-    # be optimized for the different survey programs
-    scale_x = scale
-    scale_y = scale
-    # set the pixelsscale to 0.1 arcsec/pixel
-    pixelscale = models.Scale(scale_x / 3600.0) & models.Scale(scale_y / 3600.0)
+    bounding_box : None or 4-tuple
+        The bounding box in detector/pixel space. Form of input is:
+        (x_left, x_right, y_bottom, y_top)
 
-    pixelshift = models.Shift(-1.0 * shiftx) & models.Shift(-1.0 * shifty)
+    name : str
+        Value of the `name` attribute of the GWCS object.
+
+    Returns
+    -------
+    wcs : wcs.GWCS
+        The GWCS object created.
+    """
+    pixelshift = models.Shift(-wcsinfo["x_ref"], name="crpix1") & models.Shift(
+        -wcsinfo["y_ref"], name="crpix2"
+    )
+    pixelscale = models.Scale(wcsinfo["pixel_scale"], name="cdelt1") & models.Scale(
+        wcsinfo["pixel_scale"], name="cdelt2"
+    )
     tangent_projection = models.Pix2Sky_TAN()
-    celestial_rotation = models.RotateNative2Celestial(ra_center, dec_center, 180.0)
-    det2sky = pixelshift | pixelscale | tangent_projection | celestial_rotation
+    celestial_rotation = models.RotateNative2Celestial(
+        wcsinfo["ra_ref"], wcsinfo["dec_ref"], 180.0
+    )
+
+    matrix = wcsinfo.get("rotation_matrix", None)
+    if matrix:
+        matrix = np.array(matrix)
+    else:
+        orientat = wcsinfo.get("orientat", 0.0)
+        matrix = wcs_util.calc_rotation_matrix(
+            np.deg2rad(orientat), v3i_yangle=0.0, vparity=1
+        )
+        matrix = np.reshape(matrix, (2, 2))
+    rotation = models.AffineTransformation2D(matrix, name="pc_rotation_matrix")
+    det2sky = (
+        pixelshift | rotation | pixelscale | tangent_projection | celestial_rotation
+    )
+
     detector_frame = coordinate_frames.Frame2D(
         name="detector", axes_names=("x", "y"), unit=(u.pix, u.pix)
     )
     sky_frame = coordinate_frames.CelestialFrame(
         reference_frame=coordinates.ICRS(), name="icrs", unit=(u.deg, u.deg)
     )
-    wcsobj = WCS([(detector_frame, det2sky), (sky_frame, None)])
-    wcsobj.bounding_box = bounding_box
+    wcsobj = WCS([(detector_frame, det2sky), (sky_frame, None)], name=name)
+
+    if bounding_box:
+        wcsobj.bounding_box = bounding_box
 
     return wcsobj