RCAL-112: Update photom step for WFI data

CHANGES.rst

@@ -12,6 +12,12 @@ Documentation
 - Add documentation for DNS build 0.5, e.g. reference array trimming [#457]
 
 
+photom
+------
+
+ - Added photom correction step and unit tests. [#469]
+
+
 0.6.0 (2022-03-02)
 ==================
 

romancal/lib/suffix.py

@@ -67,6 +67,7 @@
     'assign_wcs',
     'linearity',
     'rampfitstep',
+    'photomstep',
     'pipeline',
     'dq_init',
     'linearitystep',

romancal/photom/__init__.py

@@ -0,0 +1,3 @@
+from .photom_step import PhotomStep
+
+__all__ = ['PhotomStep']

romancal/photom/photom.py

@@ -0,0 +1,110 @@
+import logging
+import warnings
+from astropy import units as u
+
+log = logging.getLogger(__name__)
+log.setLevel(logging.DEBUG)
+
+
+def photom_io(input_model, photom_metadata):
+    """
+    Combine photometric scalar conversion factors and add to the science metadata.
+
+    Parameters
+    ----------
+    input_model : Roman level 2 image datamodel
+        Input Roman datamodel
+
+    photom_metadata : dict
+        Set of matched photom meta data keys
+
+    Returns
+    -------
+
+    """
+    # Get the scalar conversion factor.
+    conversion = photom_metadata['photmjsr']  # unit is MJy / sr
+
+    # Store the conversion factor in the meta data
+    log.info(f'photmjsr value: {conversion:.6g}')
+    input_model.meta.photometry.conversion_megajanskys = conversion
+    input_model.meta.photometry.conversion_microjanskys = conversion.to(
+        u.microjansky / u.arcsecond ** 2)
+
+    # Get the scalar conversion uncertainty factor
+    uncertainty_conv = photom_metadata['uncertainty']
+
+    # Store the uncertainty conversion factor in the meta data
+    log.info(f'uncertainty value: {conversion:.6g}')
+    input_model.meta.photometry.conversion_megajanskys_uncertainty = uncertainty_conv
+    input_model.meta.photometry.conversion_microjanskys_uncertainty = uncertainty_conv.to(
+        u.microjansky / u.arcsecond ** 2)
+
+    # Return updated input model
+    return input_model
+
+
+def save_area_info(input_model, photom_parameters):
+    """
+    Read the pixel area value in the photom parameters, then convert and
+    copy them to the metadata of the input datamodel.
+
+    Parameters
+    ----------
+    input_model : Roman level 2 image datamodel
+        Input Roman datamodel
+
+    photom_parameters : dict
+        Photom parameter object
+    """
+
+    # Load the average pixel area values from the photom reference file header
+    area_ster = photom_parameters["pixelareasr"]
+    area_a2 = photom_parameters["pixelareasr"].to(u.arcsecond**2)
+
+    # Copy the pixel area values to the input model
+    log.debug(f'pixelarea_steradians = {area_ster}, pixelarea_arcsecsq = {area_a2}')
+    input_model.meta.photometry.pixelarea_arcsecsq = area_a2
+    input_model.meta.photometry.pixelarea_steradians = area_ster
+
+    # Return updated input model
+    return input_model
+
+
+def apply_photom(input_model, photom):
+    """
+    Retrieve the conversion factor from the photom reference datamodel
+    that is appropriate to the instrument mode. The scalar factor is written
+    to the photmjsr and uncertainty keywords in the model.
+
+    For WFI, matching is based on optical_element.
+
+    Parameters
+    ----------
+    input_model : Roman level 2 image datamodel
+        Input Roman datamodel
+
+    photom : Roman Photom reference datamodel
+        Photom Roman datamodel
+
+    Returns
+    -------
+    output_model : Roman level 2 image datamodel
+        Output Roman datamodel with the flux calibration keywords updated
+    """
+    # Obtain photom parameters for the selected optical element
+    try:
+        photom_parameters = photom.phot_table[input_model.meta.instrument.optical_element.upper()]
+    except KeyError:
+        warnings.warn(f'No matching photom parameters for '
+                      f'{input_model.meta.instrument.optical_element}')
+        return input_model
+
+    # Copy pixel area information to output datamodel
+    output_model = save_area_info(input_model, photom_parameters)
+
+    # Copy conversions to output model
+    output_model = photom_io(output_model, photom_parameters)
+
+    # Return updated output model
+    return output_model

romancal/photom/photom_step.py

@@ -0,0 +1,62 @@
+#! /usr/bin/env python
+
+from romancal.stpipe import RomanStep
+from romancal.photom import photom
+import roman_datamodels as rdm
+
+__all__ = ["PhotomStep"]
+
+
+class PhotomStep(RomanStep):
+    """
+    PhotomStep: Module for loading photometric conversion information from
+        reference files and attaching to the input science data model
+    """
+
+    reference_file_types = ['photom']
+
+    def process(self, input):
+        """Perform the photometric calibration step
+
+        Parameters
+        ----------
+        input : Roman level 2 image datamodel (wfi_image-1.x.x)
+            input roman datamodel
+
+        Returns
+        -------
+        output_model : Roman level 2 image datamodel (wfi_image-1.x.x)
+            output roman datamodel
+        """
+
+        # Open the input data model
+        with rdm.open(input) as input_model:
+
+            # Get reference file
+            reffile = self.get_reference_file(input_model, "photom")
+
+            # Open the relevant photom reference file as a datamodel
+            if reffile is not None:
+                # If there is a reference file, perform photom application
+                photom_model = rdm.open(reffile)
+                self.log.debug(f'Using PHOTOM ref file: {reffile}')
+
+                # Do the correction
+                output_model = photom.apply_photom(input_model, photom_model)
+                output_model.meta.cal_step.photom = 'COMPLETE'
+
+            else:
+                # Skip Photom step if no photom file
+                self.log.warning('No PHOTOM reference file found')
+                self.log.warning('Photom step will be skipped')
+                input_model.meta.cal_step.photom = 'SKIPPED'
+                return input_model
+
+        # Close the input and reference files
+        input_model.close()
+        try:
+            photom_model.close()
+        except AttributeError:
+            pass
+
+        return output_model