Rcal 504 develop the resample step (spacetelescope#787)

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Nadia Dencheva <[email protected]>
Co-authored-by: D Davis <[email protected]>
Co-authored-by: William Jamieson <[email protected]>
Co-authored-by: Jonathan Eisenhamer <[email protected]>
Co-authored-by: Jonathan Eisenhamer <[email protected]>
Co-authored-by: Brett M. Morris <[email protected]>
Co-authored-by: Zach Burnett <[email protected]>
Co-authored-by: Zach Burnett <[email protected]>

CHANGES.rst

@@ -47,6 +47,10 @@ refpix
 
 - Update cal_step, add suffix and add to the exposure pipeline [#890]
 
+resample
+--------
+- Implement resampling step. [#787]
+
 
 0.12.0 (2023-08-18)
 ===================

docs/conf.py

@@ -67,6 +67,8 @@ def check_sphinx_version(expected_version):
     "numpy": ("https://numpy.org/devdocs", None),
     "scipy": ("http://scipy.github.io/devdocs", None),
     "matplotlib": ("http://matplotlib.org/", None),
+    "gwcs": ("https://gwcs.readthedocs.io/en/latest/", None),
+    "astropy": ("https://docs.astropy.org/en/stable/", None),
 }
 
 if sys.version_info[0] == 2:

docs/roman/package_index.rst

@@ -16,3 +16,4 @@
    tweakreg/index.rst
    outlier_detection/index.rst
    skymatch/index.rst
+   resample/index.rst

docs/roman/pipeline_steps.rst

@@ -14,3 +14,4 @@ Pipeline Steps
    photom/index.rst
    source_detection/index.rst
    tweakreg/index.rst
+   resample/index.rst

docs/roman/resample/arguments.rst

@@ -0,0 +1,104 @@
+.. _resample_step_args:
+
+Step Arguments
+==============
+The ``resample`` step has the following optional arguments that control
+the behavior of the processing and the characteristics of the resampled
+image.
+
+``--pixfrac`` (float, default=1.0)
+    The fraction by which input pixels are "shrunk" before being drizzled
+    onto the output image grid, given as a real number between 0 and 1.
+
+``--kernel`` (str, default='square')
+    The form of the kernel function used to distribute flux onto the output
+    image.  Available kernels are `square`, `gaussian`, `point`, `tophat`,
+    `turbo`, `lanczos2`, and `lanczos3`.
+
+``--pixel_scale_ratio`` (float, default=1.0)
+    Ratio of input to output pixel scale.  A value of 0.5 means the output
+    image would have 4 pixels sampling each input pixel.
+    Ignored when ``pixel_scale`` or ``output_wcs`` are provided.
+
+``--pixel_scale`` (float, default=None)
+    Absolute pixel scale in ``arcsec``. When provided, overrides
+    ``pixel_scale_ratio``. Ignored when ``output_wcs`` is provided.
+
+``--rotation`` (float, default=None)
+    Position angle of output image’s Y-axis relative to North.
+    A value of 0.0 would orient the final output image to be North up.
+    The default of `None` specifies that the images will not be rotated,
+    but will instead be resampled in the default orientation for the camera
+    with the x and y axes of the resampled image corresponding
+    approximately to the detector axes. Ignored when ``pixel_scale``
+    or ``output_wcs`` are provided.
+
+``--crpix`` (tuple of float, default=None)
+    Position of the reference pixel in the image array in the ``x, y`` order.
+    If ``crpix`` is not specified, it will be set to the center of the bounding
+    box of the returned WCS object. When supplied from command line, it should
+    be a comma-separated list of floats. Ignored when ``output_wcs``
+    is provided.
+
+``--crval`` (tuple of float, default=None)
+    Right ascension and declination of the reference pixel. Automatically
+    computed if not provided. When supplied from command line, it should be a
+    comma-separated list of floats. Ignored when ``output_wcs`` is provided.
+
+``--output_shape`` (tuple of int, default=None)
+    Shape of the image (data array) using "standard" ``nx`` first and ``ny``
+    second (as opposite to the ``numpy.ndarray`` convention - ``ny`` first and
+    ``nx`` second). This value will be assigned to
+    ``pixel_shape`` and ``array_shape`` properties of the returned
+    WCS object. When supplied from command line, it should be a comma-separated
+    list of integers ``nx, ny``.
+
+    .. note::
+        Specifying ``output_shape`` *is required* when the WCS in
+        ``output_wcs`` does not have ``bounding_box`` property set.
+
+``--output_wcs`` (str, default='')
+    File name of a ``ASDF`` file with a GWCS stored under the ``"wcs"`` key
+    under the root of the file. The output image size is determined from the
+    bounding box of the WCS (if any). Argument ``output_shape`` overrides
+    computed image size and it is required when output WCS does not have
+    ``bounding_box`` property set.
+
+    .. note::
+        When ``output_wcs`` is specified, WCS-related arguments such as
+        ``pixel_scale_ratio``, ``pixel_scale``, ``rotation``, ``crpix``,
+        and ``crval`` will be ignored.
+
+``--fillval`` (str, default='INDEF')
+    The value to assign to output pixels that have zero weight or do not
+    receive any flux from any input pixels during drizzling.
+
+``--weight_type`` (str, default='ivm')
+    The weighting type for each input image.
+    If `weight_type=ivm` (the default), the scaling value
+    will be determined per-pixel using the inverse of the read noise
+    (VAR_RNOISE) array stored in each input image. If the VAR_RNOISE array does
+    not exist, the variance is set to 1 for all pixels (equal weighting).
+    If `weight_type=exptime`, the scaling value will be set equal to the
+    exposure time found in the image header.
+
+``--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/resample/index.rst

@@ -0,0 +1,16 @@
+.. _resample_step:
+
+==========
+Resample
+==========
+
+.. toctree::
+   :maxdepth: 2
+
+   main.rst
+   arguments.rst
+   resample_step.rst
+   resample.rst
+   resample_utils.rst
+
+.. automodapi:: romancal.resample

docs/roman/resample/main.rst

@@ -0,0 +1,118 @@
+Description
+===========
+
+:Classes: `romancal.resample.ResampleStep`
+:Alias: resample
+
+This routine will resample each input 2D image based on the WCS and
+distortion information, and will combine multiple resampled images
+into a single undistorted product.  The distortion information should have
+been incorporated into the image using the :ref:`assign_wcs <assign_wcs_step>`
+step.
+
+The ``resample`` step can take:
+
+  * a single 2D input image (in the format of either a string with the full
+    path and filename of an ASDF file or a Roman
+    Datamodel/:py:class:`~romancal.datamodels.container.ModelContainer`);
+  * an association table (in JSON format).
+
+The parameters for the drizzle operation itself are set by
+:py:func:`~romancal.resample.resample_step.ResampleStep.set_drizzle_defaults`.
+The exact values used depends on the number of input images being combined
+and the filter being used. Other information may be added as selection criteria
+later, but during the :py:class:`~romancal.resample.resample_step.ResampleStep`
+instantiation, only basic information is set.
+
+The output product is determined by using the WCS information of all inputs,
+even if it is just a single image. The output WCS defines a
+field-of-view that encompasses the undistorted footprints on the sky
+of all the input images with the same orientation and plate scale
+as the first listed input image.
+
+This step uses the interface to the C-based `cdriz` routine to do the
+resampling via the
+drizzle method (`Fruchter and Hook, PASP 2002`_).
+The input-to-output pixel mapping is determined via a mapping function
+derived from the WCS of each input image and the WCS of the defined
+output product. The mapping function is created by
+:py:func:`~romancal.resample.resample_utils.reproject` and passed on to
+`cdriz` to drive the actual drizzling to create the output product.
+
+Context Image
+-------------
+
+In addition to the resampled image data, resample step also creates a
+"context image" stored in the ``con`` attribute in the output data model.
+Each pixel in the context image is a bit field that encodes
+information about which input image has contributed to the corresponding
+pixel in the resampled data array. Context image uses 32 bit integers to encode
+this information and hence it can keep track of 32 input images at most.
+
+For any given pixel, the first bit corresponds to the first input image,
+the second bit corrsponds to the second input image, and so on.
+If the number of input images is larger than 32, then it is necessary to
+have multiple context images ("planes") to hold information about all input
+images with the first plane encoding which of the first 32 images
+(indexed from 0 through 32) contributed to the output data pixel, second plane
+representing next 32 input images (indexed from 33 through 64), etc.
+For this reason, the context image is a 3D array of type `numpy.int32` and shape
+``(np, ny, nx)``, where ``nx`` and ``ny`` are the dimensions of the image's data
+and ``np`` is the number of "planes", which is equal to
+``(number of input images - 1) // 32 + 1``. If a bit at position ``k`` in a pixel
+with coordinates ``(p, y, x)`` is 0 then input image number ``32 * p + k``
+(0-indexed) did not contribute to the output data pixel with array coordinates
+``(y, x)``, and if that bit is 1 then input image number ``32 * p + k`` did
+contribute to the pixel ``(y, x)`` in the resampled image.
+
+As an example, let's assume we have 8 input images. Then, when ``con`` pixel
+values are displayed using binary representation (and decimal in parenthesis),
+one could see values like this::
+
+    00000001 (1) - only first input image contributed to this output pixel;
+    00000010 (2) - 2nd input image contributed;
+    00000100 (4) - 3rd input image contributed;
+    10000000 (128) - 8th input image contributed;
+    10000100 (132=128+4) - 3rd and 8th input images contributed;
+    11001101 (205=1+4+8+64+128) - input images 1, 3, 4, 7, 8 have contributed
+    to this output pixel.
+
+In order to test if a specific input image contributed to an output pixel,
+one needs to use bitwise operations. Using the example above, to test whether
+input images number 4 and 5 have contributed to the output pixel whose
+corresponding ``con`` value is 205 (11001101 in binary form) we can do
+the following:
+
+>>> bool(205 & (1 << (5 - 1)))  # (205 & 16) = 0 (== 0 => False): did NOT contribute
+False
+>>> bool(205 & (1 << (4 - 1)))  # (205 & 8) = 8 (!= 0 => True): did contribute
+True
+
+In general, to get a list of all input images that have contributed to an
+output resampled pixel with image coordinates ``(x, y)``, and given a
+context array ``con``, one can do something like this:
+
+.. doctest-skip::
+
+    >>> import numpy as np
+    >>> np.flatnonzero([v & (1 << k) for v in con[:, y, x] for k in range(32)])
+
+For convenience, this functionality was implemented in the
+:py:func:`~romancal.resample.resample_utils.decode_context` function.
+
+
+References
+----------
+
+* `Fruchter and Hook, PASP 2002`_: full description of the drizzling algorithm.
+
+* `Casertano et al., AJ 2000`_ (Appendix A2): description of the inverse variance
+  map method.
+
+* `DrizzlePac Handbook`_: description of the drizzle parameters and other useful
+  drizzle-related resources.
+
+
+.. _Fruchter and Hook, PASP 2002: https://doi.org/10.1086/338393
+.. _Casertano et al., AJ 2000: https://doi.org/10.1086/316851
+.. _DrizzlePac Handbook: http://drizzlepac.stsci.edu

docs/roman/resample/resample.rst

@@ -0,0 +1,6 @@
+.. resample_:
+
+Python Interface to Drizzle: ResampleData()
+===========================================
+
+.. automodapi:: romancal.resample.resample

docs/roman/resample/resample_step.rst

@@ -0,0 +1,6 @@
+.. resample_step_:
+
+Python Step Interface: ResampleStep()
+=====================================
+
+.. automodapi:: romancal.resample.resample_step

docs/roman/resample/resample_utils.rst

@@ -0,0 +1,7 @@
+Resample Utilities
+==================
+
+.. currentmodule:: romancal.resample.resample_utils
+
+.. automodule:: romancal.resample.resample_utils
+   :members:

pyproject.toml

@@ -33,6 +33,7 @@ dependencies = [
     'tweakwcs >=0.8.0',
     'spherical-geometry >= 1.2.22',
     'stsci.imagestats >= 1.6.3',
+    'drizzle >= 1.13.7',
     'webbpsf == 1.1.1',
 ]
 dynamic = ['version']
@@ -46,6 +47,7 @@ docs = [
     'sphinx-automodapi',
     'sphinx-rtd-theme',
     'stsci-rtd-theme',
+    'sphinx-autobuild',
     'tomli; python_version <="3.11"',
 ]
 test = [

romancal/assign_wcs/assign_wcs_step.py

@@ -12,7 +12,7 @@
 
 from ..stpipe import RomanStep
 from . import pointing
-from .utils import wcs_bbox_from_shape
+from .utils import add_s_region, wcs_bbox_from_shape
 
 log = logging.getLogger(__name__)
 log.setLevel(logging.DEBUG)
@@ -66,10 +66,9 @@ def load_wcs(input_model, reference_files=None):
 
     if reference_files is not None:
         for ref_type, ref_file in reference_files.items():
-            if ref_file not in ["N/A", ""]:
-                reference_files[ref_type] = ref_file
-            else:
-                reference_files[ref_type] = None
+            reference_files[ref_type] = (
+                ref_file if ref_file not in ["N/A", ""] else None
+            )
     else:
         reference_files = {}
 
@@ -87,7 +86,11 @@ def load_wcs(input_model, reference_files=None):
     distortion = wfi_distortion(output_model, reference_files)
     tel2sky = pointing.v23tosky(output_model)
 
-    pipeline = [Step(detector, distortion), Step(v2v3, tel2sky), Step(world, None)]
+    pipeline = [
+        Step(detector, distortion),
+        Step(v2v3, tel2sky),
+        Step(world, None),
+    ]
     wcs = WCS(pipeline)
     if wcs.bounding_box is None:
         wcs.bounding_box = wcs_bbox_from_shape(output_model.data.shape)
@@ -137,49 +140,3 @@ def wfi_distortion(model, reference_files):
         transform.bounding_box = bbox
 
     return transform
-
-
-def add_s_region(model):
-    """
-    Calculate the detector's footprint using ``WCS.footprint`` and save it in the
-    ``S_REGION`` keyword
-
-    Parameters
-    ----------
-    model : `~roman_datamodels.datamodels.ImageModel`
-        The data model for processing
-
-    Returns
-    -------
-    A formatted string representing the detector's footprint
-    """
-
-    bbox = model.meta.wcs.bounding_box
-
-    if bbox is None:
-        bbox = wcs_bbox_from_shape(model.data.shape)
-
-    # footprint is an array of shape (2, 4) - i.e. 4 values for RA and 4 values for
-    # Dec - as we are interested only in the footprint on the sky
-    footprint = model.meta.wcs.footprint(bbox, center=True, axis_type="spatial").T
-    # take only imaging footprint
-    footprint = footprint[:2, :]
-
-    # Make sure RA values are all positive
-    negative_ind = footprint[0] < 0
-    if negative_ind.any():
-        footprint[0][negative_ind] = 360 + footprint[0][negative_ind]
-
-    footprint = footprint.T
-    update_s_region_keyword(model, footprint)
-
-
-def update_s_region_keyword(model, footprint):
-    s_region = "POLYGON ICRS " + " ".join([str(x) for x in footprint.ravel()]) + " "
-    log.info(f"S_REGION VALUES: {s_region}")
-    if "nan" in s_region:
-        # do not update s_region if there are NaNs.
-        log.info("There are NaNs in s_region, S_REGION not updated.")
-    else:
-        model.meta.wcsinfo.s_region = s_region
-        log.info(f"Update S_REGION to {model.meta.wcsinfo.s_region}")