RCAL-911 & 932: remove units from MOS and ELP pipelines.

changes/1445.general.rst

@@ -0,0 +1 @@
+Remove units from romancal.

docs/conf.py

@@ -273,7 +273,7 @@ def check_sphinx_version(expected_version):
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the
 # documentation.
-html_theme_options = {"collapse_navigation": True, "display_version": True}
+html_theme_options = {"collapse_navigation": True}
 #        "nosidebar": "false",
 #        "sidebarbgcolor": "#4db8ff",
 #        "sidebartextcolor": "black",

pyproject.toml

@@ -24,9 +24,9 @@ dependencies = [
     "pyparsing >=2.4.7",
     "requests >=2.26",
     # "rad>=0.21.0,<0.22.0",
-    "rad @ git+https://github.com/spacetelescope/rad.git",
+    "rad @ git+https://github.com/mairanteodoro/rad.git@RCAL-932-remove-units-in-exposure-level-pipeline",
     # "roman_datamodels>=0.21.0,<0.22.0",
-    "roman_datamodels @ git+https://github.com/spacetelescope/roman_datamodels.git",
+    "roman_datamodels @ git+https://github.com/mairanteodoro/roman_datamodels.git@RCAL-932-remove-units-in-exposure-level-pipeline",
     "scipy >=1.11",
 #    "stcal>=1.8.0,<1.9.0",
     "stcal @ git+https://github.com/spacetelescope/stcal.git@main",

romancal/dark_current/dark_current_step.py

@@ -54,7 +54,7 @@ def process(self, input):
             # Do the dark correction
             out_model = input_model
             nresultants = len(input_model.meta.exposure["read_pattern"])
-            out_model.data -= dark_model.data[:nresultants]
+            out_model.data -= dark_model.data[:nresultants].value
             out_model.pixeldq |= dark_model.dq
             out_model.meta.cal_step.dark = "COMPLETE"
 

romancal/dark_current/tests/test_dark.py

@@ -5,7 +5,6 @@
 import numpy as np
 import pytest
 import roman_datamodels as rdm
-from astropy import units as u
 from roman_datamodels import maker_utils
 from roman_datamodels.datamodels import DarkRefModel, RampModel
 
@@ -60,20 +59,20 @@ def test_dark_step_subtraction(instrument, exptype):
 
     # populate data array of science cube
     for i in range(0, 20):
-        ramp_model.data[0, 0, i] = i * ramp_model.data.unit
+        ramp_model.data[0, 0, i] = i
         darkref_model.data[0, 0, i] = i * 0.1 * darkref_model.data.unit
     orig_model = ramp_model.copy()
 
     # Perform Dark Current subtraction step
     result = DarkCurrentStep.call(ramp_model, override_dark=darkref_model)
 
     # check that the dark file is subtracted frame by frame from the science data
-    diff = orig_model.data.value - darkref_model.data.value
+    diff = orig_model.data - darkref_model.data.value
 
     # test that the output data file is equal to the difference found when subtracting
     # reffile from sci file
     np.testing.assert_array_equal(
-        result.data.value, diff, err_msg="dark file should be subtracted from sci file "
+        result.data, diff, err_msg="dark file should be subtracted from sci file "
     )
 
 
@@ -148,7 +147,7 @@ def create_ramp_and_dark(shape, instrument, exptype):
     ramp.meta.instrument.detector = "WFI01"
     ramp.meta.instrument.optical_element = "F158"
     ramp.meta.exposure.type = exptype
-    ramp.data = u.Quantity(np.ones(shape, dtype=np.float32), u.DN, dtype=np.float32)
+    ramp.data = np.ones(shape, dtype=np.float32)
     ramp_model = RampModel(ramp)
 
     # Create dark model

romancal/dq_init/tests/test_dq_init.py

@@ -1,6 +1,5 @@
 import numpy as np
 import pytest
-from astropy import units as u
 from roman_datamodels import maker_utils, stnode
 from roman_datamodels.datamodels import MaskRefModel, ScienceRawModel
 from roman_datamodels.dqflags import pixel
@@ -199,9 +198,7 @@ def test_dqinit_step_interface(instrument, exptype):
     wfi_sci_raw.meta["guidestar"]["gw_window_xstart"] = 1012
     wfi_sci_raw.meta["guidestar"]["gw_window_xsize"] = 16
     wfi_sci_raw.meta.exposure.type = exptype
-    wfi_sci_raw.data = u.Quantity(
-        np.ones(shape, dtype=np.uint16), u.DN, dtype=np.uint16
-    )
+    wfi_sci_raw.data = np.ones(shape, dtype=np.uint16)
     wfi_sci_raw[extra_key] = extra_value
     wfi_sci_raw_model = ScienceRawModel(wfi_sci_raw)
 
@@ -251,9 +248,7 @@ def test_dqinit_refpix(instrument, exptype):
     wfi_sci_raw.meta["guidestar"]["gw_window_xstart"] = 1012
     wfi_sci_raw.meta["guidestar"]["gw_window_xsize"] = 16
     wfi_sci_raw.meta.exposure.type = exptype
-    wfi_sci_raw.data = u.Quantity(
-        np.ones(shape, dtype=np.uint16), u.DN, dtype=np.uint16
-    )
+    wfi_sci_raw.data = np.ones(shape, dtype=np.uint16)
     wfi_sci_raw_model = ScienceRawModel(wfi_sci_raw)
 
     # Create mask model
@@ -272,7 +267,7 @@ def test_dqinit_refpix(instrument, exptype):
 
     # check if reference pixels are correct
     assert result.data.shape == (2, 20, 20)  # no pixels should be trimmed
-    assert result.amp33.value.shape == (2, 4096, 128)
+    assert result.amp33.shape == (2, 4096, 128)
     assert result.border_ref_pix_right.shape == (2, 20, 4)
     assert result.border_ref_pix_left.shape == (2, 20, 4)
     assert result.border_ref_pix_top.shape == (2, 4, 20)
@@ -304,9 +299,7 @@ def test_dqinit_resultantdq(instrument, exptype):
     wfi_sci_raw.meta["guidestar"]["gw_window_xsize"] = 16
     wfi_sci_raw.meta.exposure.type = exptype
     wfi_sci_raw.resultantdq[1, 12, 12] = pixel["DROPOUT"]
-    wfi_sci_raw.data = u.Quantity(
-        np.ones(shape, dtype=np.uint16), u.DN, dtype=np.uint16
-    )
+    wfi_sci_raw.data = np.ones(shape, dtype=np.uint16)
     wfi_sci_raw_model = ScienceRawModel(wfi_sci_raw)
 
     # Create mask model
@@ -354,9 +347,7 @@ def test_dqinit_getbestref(instrument, exptype):
     wfi_sci_raw.meta["guidestar"]["gw_window_xstart"] = 1012
     wfi_sci_raw.meta["guidestar"]["gw_window_xsize"] = 16
     wfi_sci_raw.meta.exposure.type = exptype
-    wfi_sci_raw.data = u.Quantity(
-        np.ones(shape, dtype=np.uint16), u.DN, dtype=np.uint16
-    )
+    wfi_sci_raw.data = np.ones(shape, dtype=np.uint16)
     wfi_sci_raw_model = ScienceRawModel(wfi_sci_raw)
 
     # Perform Data Quality application step

romancal/flatfield/flat_field.py

@@ -5,7 +5,6 @@
 import logging
 
 import numpy as np
-from astropy import units as u
 from roman_datamodels.dqflags import pixel
 
 log = logging.getLogger(__name__)
@@ -105,9 +104,7 @@ def apply_flat_field(science, flat):
     flat_data[np.where(flat_bad)] = 1.0
     # Now let's apply the correction to science data and error arrays.  Rely
     # on array broadcasting to handle the cubes
-    science.data = u.Quantity(
-        (science.data.value / flat_data), u.DN / u.s, dtype=science.data.dtype
-    )
+    science.data = (science.data / flat_data).astype(science.data.dtype)
 
     # Update the variances using BASELINE algorithm.  For guider data, it has
     # not gone through ramp fitting so there is no Poisson noise or readnoise
@@ -121,8 +118,6 @@ def apply_flat_field(science, flat):
         science.var_flat = science.data**2 / flat_data_squared * flat_err**2
 
     science.err = np.sqrt(science.var_poisson + science.var_rnoise + science.var_flat)
-    science.err = science.err.to(science.data.unit)
-    # Workaround for https://github.com/astropy/astropy/issues/16055
 
     # Combine the science and flat DQ arrays
     science.dq = np.bitwise_or(science.dq, flat_dq)

romancal/flatfield/tests/test_flatfield.py

@@ -1,6 +1,5 @@
 import numpy as np
 import pytest
-from astropy import units as u
 from astropy.time import Time
 from roman_datamodels import maker_utils, stnode
 from roman_datamodels.datamodels import FlatRefModel, ImageModel
@@ -26,22 +25,12 @@ def test_flatfield_step_interface(instrument, exptype):
     wfi_image.meta.instrument.detector = "WFI01"
     wfi_image.meta.instrument.optical_element = "F158"
     wfi_image.meta.exposure.type = exptype
-    wfi_image.data = u.Quantity(
-        np.ones(shape, dtype=np.float32), u.DN / u.s, dtype=np.float32
-    )
+    wfi_image.data = np.ones(shape, dtype=np.float32)
     wfi_image.dq = np.zeros(shape, dtype=np.uint32)
-    wfi_image.err = u.Quantity(
-        np.zeros(shape, dtype=np.float32), u.DN / u.s, dtype=np.float32
-    )
-    wfi_image.var_poisson = u.Quantity(
-        np.zeros(shape, dtype=np.float32), u.DN**2 / u.s**2, dtype=np.float32
-    )
-    wfi_image.var_rnoise = u.Quantity(
-        np.zeros(shape, dtype=np.float32), u.DN**2 / u.s**2, dtype=np.float32
-    )
-    wfi_image.var_flat = u.Quantity(
-        np.zeros(shape, dtype=np.float32), u.DN**2 / u.s**2, dtype=np.float32
-    )
+    wfi_image.err = np.zeros(shape, dtype=np.float32)
+    wfi_image.var_poisson = np.zeros(shape, dtype=np.float32)
+    wfi_image.var_rnoise = np.zeros(shape, dtype=np.float32)
+    wfi_image.var_flat = np.zeros(shape, dtype=np.float32)
 
     wfi_image_model = ImageModel(wfi_image)
     flatref = stnode.FlatRef()

romancal/flux/flux_step.py

@@ -2,7 +2,6 @@
 
 import logging
 
-import astropy.units as u
 from roman_datamodels import datamodels
 
 from ..datamodels import ModelLibrary
@@ -14,10 +13,6 @@
 __all__ = ["FluxStep"]
 
 
-# Define expected Level 2 units
-LV2_UNITS = u.DN / u.s
-
-
 class FluxStep(RomanStep):
     """Apply flux scaling to count-rate data
 
@@ -104,26 +99,19 @@
     DATA = ("data", "err")
     VARIANCES = ("var_rnoise", "var_poisson", "var_flat")
 
-    if model.data.unit == model.meta.photometry.conversion_megajanskys.unit:
-        log.info(
-            f"Input data is already in flux units of {model.meta.photometry.conversion_megajanskys.unit}."
-        )
-        log.info("Flux correction already applied.")
-        return
-
-    if model.data.unit != LV2_UNITS:
+    if model.meta.cal_step["flux"] == "COMPLETE":
         message = (
-            f"Input data units {model.data.unit} are not in the expected units of {LV2_UNITS}"
-            "\nAborting flux correction"
+            "Input data is already in flux units of MJy/sr."
+            "\nFlux correction already applied."
         )
-        [log.error(line) for line in message.splitlines()]
-        raise ValueError(message)
+        log.info(message)
+        return
 
     # Apply the correction.
     # The end goal in units is to have MJy/sr. The scale is in MJy/sr also.
     # Hence the extra factor of s/DN must be applied to cancel DN/s.
     log.debug("Flux correction being applied")
-    c_mj = model.meta.photometry.conversion_megajanskys / model.data.unit
+    c_mj = model.meta.photometry.conversion_megajanskys
     for data in DATA:
         model[data] = model[data] * c_mj
     for variance in VARIANCES:

romancal/flux/tests/test_flux_step.py

@@ -7,7 +7,6 @@
 
 from romancal.datamodels import ModelLibrary
 from romancal.flux import FluxStep
-from romancal.flux.flux_step import LV2_UNITS
 
 
 @pytest.mark.parametrize(
@@ -23,7 +22,7 @@
 def test_attributes(flux_step, attr, factor):
     """Test that the attribute has been scaled by the right factor"""
     original, result = flux_step
-    c_unit = 1.0 / LV2_UNITS
+    c_unit = 1.0
 
     # Handle difference between just a single image and a list.
     if isinstance(original, datamodels.ImageModel):
@@ -86,31 +85,15 @@ def flux_step(request):
 @pytest.fixture(scope="module")
 def image_model():
     """Product a basic ImageModel"""
-    # Create a random image and specify a conversion.
+    # Create a random image and specify a conversion
     rng = np.random.default_rng()
     shape = (10, 10)
     image_model = maker_utils.mk_datamodel(datamodels.ImageModel, shape=shape)
-    image_model.data = u.Quantity(
-        rng.poisson(2.5, size=shape).astype(np.float32),
-        LV2_UNITS,
-        dtype=np.float32,
-    )
-    image_model.var_rnoise = u.Quantity(
-        rng.normal(1, 0.05, size=shape).astype(np.float32),
-        LV2_UNITS**2,
-        dtype=np.float32,
-    )
-    image_model.var_poisson = u.Quantity(
-        rng.poisson(1, size=shape).astype(np.float32),
-        LV2_UNITS**2,
-        dtype=np.float32,
-    )
-    image_model.var_flat = u.Quantity(
-        rng.uniform(0, 1, size=shape).astype(np.float32),
-        LV2_UNITS**2,
-        dtype=np.float32,
-    )
-    image_model.meta.photometry.conversion_megajanskys = 2.0 * u.MJy / u.sr
+    image_model.data = rng.poisson(2.5, size=shape).astype(np.float32)
+    image_model.var_rnoise = rng.normal(1, 0.05, size=shape).astype(np.float32)
+    image_model.var_poisson = rng.poisson(1, size=shape).astype(np.float32)
+    image_model.var_flat = rng.uniform(0, 1, size=shape).astype(np.float32)
+    image_model.meta.photometry.conversion_megajanskys = (2.0 * u.MJy / u.sr).value
 
     return image_model
 
@@ -129,6 +112,6 @@ def input_modellibrary(image_model):
     # Create and return a ModelLibrary
     image_model1 = image_model.copy()
     image_model2 = image_model.copy()
-    image_model2.meta.photometry.conversion_megajanskys = 0.5 * u.MJy / u.sr
+    image_model2.meta.photometry.conversion_megajanskys = (0.5 * u.MJy / u.sr).value
     container = ModelLibrary([image_model1, image_model2])
     return container

romancal/jump/jump_step.py

@@ -53,10 +53,10 @@ def process(self, input):
 
         # Extract the needed info from the Roman Data Model
         meta = input_model.meta
-        r_data = input_model.data.value
+        r_data = input_model.data
         r_gdq = input_model.groupdq
         r_pdq = input_model.pixeldq
-        r_err = input_model.err.value
+        r_err = input_model.err
         result = input_model
 
         # If the ramp fitting jump detection is enabled, then skip this step
@@ -106,6 +106,7 @@ def process(self, input):
             self.log.info("Maximum cores to use = %s", max_cores)
 
         # Get the gain and readnoise reference files
+        # TODO: remove units from gain and RN reference files
         gain_filename = self.get_reference_file(input_model, "gain")
         self.log.info("Using GAIN reference file: %s", gain_filename)
         gain_model = rdd.GainRefModel(gain_filename)

romancal/jump/tests/test_jump_step.py

@@ -106,10 +106,10 @@ def _setup(
         dm_ramp.meta.instrument.name = "WFI"
         dm_ramp.meta.instrument.optical_element = "F158"
 
-        dm_ramp.data = u.Quantity(data + 6.0, u.DN, dtype=np.float32)
+        dm_ramp.data = data + 6.0
         dm_ramp.pixeldq = pixdq
         dm_ramp.groupdq = gdq
-        dm_ramp.err = u.Quantity(err, u.DN, dtype=np.float32)
+        dm_ramp.err = err
 
         dm_ramp.meta.exposure.type = "WFI_IMAGE"
         dm_ramp.meta.exposure.group_time = deltatime
@@ -148,7 +148,7 @@ def test_one_CR(generate_wfi_reffiles, max_cores, setup_inputs):
     )
 
     for i in range(ngroups):
-        model1.data[i, :, :] = deltaDN * i * model1.data.unit
+        model1.data[i, :, :] = deltaDN * i
 
     first_CR_group_locs = [x for x in range(1, 7) if x % 5 == 0]
 
@@ -161,8 +161,7 @@ def test_one_CR(generate_wfi_reffiles, max_cores, setup_inputs):
     for i in range(len(CR_x_locs)):
         CR_group = next(CR_pool)
         model1.data[CR_group:, CR_y_locs[i], CR_x_locs[i]] = (
-            model1.data[CR_group:, CR_y_locs[i], CR_x_locs[i]]
-            + 5000.0 * model1.data.unit
+            model1.data[CR_group:, CR_y_locs[i], CR_x_locs[i]] + 5000.0
         )
 
     out_model = JumpStep.call(
@@ -203,7 +202,7 @@ def test_two_CRs(generate_wfi_reffiles, max_cores, setup_inputs):
     )
 
     for i in range(ngroups):
-        model1.data[i, :, :] = deltaDN * i * model1.data.unit
+        model1.data[i, :, :] = deltaDN * i
 
     first_CR_group_locs = [x for x in range(1, 7) if x % 5 == 0]
     CR_locs = [x for x in range(xsize * ysize) if x % CR_fraction == 0]
@@ -215,11 +214,10 @@ def test_two_CRs(generate_wfi_reffiles, max_cores, setup_inputs):
         CR_group = next(CR_pool)
 
         model1.data[CR_group:, CR_y_locs[i], CR_x_locs[i]] = (
-            model1.data[CR_group:, CR_y_locs[i], CR_x_locs[i]] + 5000 * model1.data.unit
+            model1.data[CR_group:, CR_y_locs[i], CR_x_locs[i]] + 5000
         )
         model1.data[CR_group + 8 :, CR_y_locs[i], CR_x_locs[i]] = (
-            model1.data[CR_group + 8 :, CR_y_locs[i], CR_x_locs[i]]
-            + 700 * model1.data.unit
+            model1.data[CR_group + 8 :, CR_y_locs[i], CR_x_locs[i]] + 700
         )
 
     out_model = JumpStep.call(