Generate L3 exposure times from resampled exposure time images

CHANGES.rst

@@ -58,6 +58,8 @@ resample
 --------
 - Implement resampling step. [#787]
 
+- Use resampled exposure time images to compute image exposure times.  [#959]
+
 stpipe
 ------
 

romancal/resample/resample.py

@@ -155,7 +155,15 @@
         )
 
         # update meta data and wcs
-        self.blank_output.meta = dict(input_models[0].meta._data.items())
+        input_model_0 = input_models[0]
+        # note we have made this input_model_0 variable so that if
+        # meta includes lazily-loaded objects, that we can successfully
+        # copy them into the metadata.  Directly running input_models[0].meta
+        # below can lead to input_models[0] going out of scope after
+        # meta is loaded but before the dictionary is constructed,
+        # which can lead to seek on closed file errors if
+        # meta contains lazily loaded objects.
+        self.blank_output.meta = dict(input_model_0.meta._data.items())
         self.blank_output.meta.wcs = self.output_wcs
 
         self.output_models = ModelContainer()
@@ -298,6 +306,9 @@
         self.resample_variance_array("var_poisson", output_model)
         self.resample_variance_array("var_flat", output_model)
 
+        # Make exposure time image
+        exptime_tot = self.resample_exposure_time(output_model)
+
         # TODO: fix unit here
         output_model.err = u.Quantity(
             np.sqrt(
@@ -313,7 +324,7 @@
             unit=output_model.err.unit,
         )
 
-        self.update_exposure_times(output_model)
+        self.update_exposure_times(output_model, exptime_tot)
 
         # TODO: fix RAD to expect a context image datatype of int32
         output_model.context = output_model.context.astype(np.uint32)
@@ -401,20 +412,82 @@
 
         setattr(output_model, name, output_variance)
 
-    def update_exposure_times(self, output_model):
+    def resample_exposure_time(self, output_model):
+        """Resample the exposure time from ``self.input_models`` to the ``output_model``.
+
+        Create an exposure time image that is the drizzled sum of the input
+        images.
+        """
+        output_wcs = output_model.meta.wcs
+        exptime_tot = np.zeros(output_model.data.shape, dtype="f4")
+
+        log.info("Resampling exposure time")
+        for model in self.input_models:
+            exptime = np.full(
+                model.data.shape, model.meta.exposure.effective_exposure_time
+            )
+
+            # create a unit weight map for all the input pixels with science data
+            inwht = resample_utils.build_driz_weight(
+                model, weight_type=None, good_bits=self.good_bits
+            )
+
+            resampled_exptime = np.zeros_like(output_model.data)
+            outwht = np.zeros_like(output_model.data)
+            outcon = np.zeros_like(output_model.context, dtype="i4")
+            # drizzle wants an i4, but datamodels wants a u4.
+
+            xmin, xmax, ymin, ymax = resample_utils.resample_range(
+                exptime.shape, model.meta.wcs.bounding_box
+            )
+
+            # resample the exptime array
+            self.drizzle_arrays(
+                exptime * u.s,  # drizzle_arrays expects these to have units
+                inwht,
+                model.meta.wcs,
+                output_wcs,
+                resampled_exptime,
+                outwht,
+                outcon,
+                pixfrac=1,  # for exposure time images, always use pixfrac = 1
+                kernel=self.kernel,
+                fillval=0,
+                xmin=xmin,
+                xmax=xmax,
+                ymin=ymin,
+                ymax=ymax,
+            )
+
+            exptime_tot += resampled_exptime.value
+
+        return exptime_tot
+
+    def update_exposure_times(self, output_model, exptime_tot):
         """Update exposure time metadata (in-place)."""
-        total_exposure_time = 0.0
+        m = exptime_tot > 0
+        if np.any(m):
+            total_exposure_time = np.median(exptime_tot[m])
+        else:
+            total_exposure_time = 0
+        max_exposure_time = np.max(exptime_tot)
+        log.info(
+            f"Mean, max exposure times: {total_exposure_time:.1f}, "
+            f"{max_exposure_time:.1f}"
+        )
         exposure_times = {"start": [], "end": []}
         for exposure in self.input_models.models_grouped:
-            total_exposure_time += exposure[0].meta.exposure.exposure_time
             exposure_times["start"].append(exposure[0].meta.exposure.start_time)
             exposure_times["end"].append(exposure[0].meta.exposure.end_time)
 
         # Update some basic exposure time values based on output_model
         output_model.meta.exposure.exposure_time = total_exposure_time
         output_model.meta.exposure.start_time = min(exposure_times["start"])
         output_model.meta.exposure.end_time = max(exposure_times["end"])
-        output_model.meta.resample.product_exposure_time = total_exposure_time
+        output_model.meta.resample.product_exposure_time = max_exposure_time
+        # we haven't filled out the L3 data model enough to put max_exposure_time
+        # somewhere sensible; I'm just dumping it in product_exposure_time
+        # for the moment.
 
     @staticmethod
     def drizzle_arrays(

romancal/resample/tests/test_resample.py

@@ -36,7 +36,10 @@ def create_image(self):
             **{
                 "meta": {
                     "wcsinfo": {"ra_ref": 10, "dec_ref": 0, "vparity": -1},
-                    "exposure": {"exposure_time": 152.04000000000002},
+                    "exposure": {
+                        "exposure_time": 152.04000000000002,
+                        "effective_exposure_time": 3.04 * 6 * 8,
+                    },
                     "observation": {
                         "program": "00005",
                         "execution_plan": 1,
@@ -458,19 +461,25 @@ def test_update_exposure_times_different_sca_same_exposure(exposure_1):
     input_models = ModelContainer(exposure_1)
     resample_data = ResampleData(input_models)
 
-    output_model = resample_data.blank_output.copy()
-    output_model.meta["resample"] = maker_utils.mk_resample()
+    output_model = resample_data.resample_many_to_one()[0]
 
-    resample_data.update_exposure_times(output_model)
+    exptime_tot = resample_data.resample_exposure_time(output_model)
+    resample_data.update_exposure_times(output_model, exptime_tot)
 
-    assert (
+    # these three SCAs overlap, so the max exposure time is 3x.
+    # get this time within 0.1 s.
+    time_difference = (
         output_model.meta.resample.product_exposure_time
-        == exposure_1[0].meta.exposure.exposure_time
+        - 3 * exposure_1[0].meta.exposure.effective_exposure_time
     )
-    assert (
+    assert np.abs(time_difference) < 0.1
+    # the median ends up being 2 exposures
+    # get this time within 0.1 s.
+    time_difference = (
         output_model.meta.exposure.exposure_time
-        == exposure_1[0].meta.exposure.exposure_time
+        - 2 * exposure_1[0].meta.exposure.effective_exposure_time
     )
+    assert np.abs(time_difference) < 0.1
     assert (
         output_model.meta.exposure.start_time == exposure_1[0].meta.exposure.start_time
     )
@@ -483,16 +492,20 @@ def test_update_exposure_times_same_sca_different_exposures(exposure_1, exposure
     input_models = ModelContainer([exposure_1[0], exposure_2[0]])
     resample_data = ResampleData(input_models)
 
-    output_model = resample_data.blank_output.copy()
-    output_model.meta["resample"] = maker_utils.mk_resample()
+    output_model = resample_data.resample_many_to_one()[0]
 
-    resample_data.update_exposure_times(output_model)
+    exptime_tot = resample_data.resample_exposure_time(output_model)
+    resample_data.update_exposure_times(output_model, exptime_tot)
 
     assert len(resample_data.input_models.models_grouped) == 2
 
-    assert output_model.meta.exposure.exposure_time == sum(
-        x.meta.exposure.exposure_time for x in input_models
+    # these exposures overlap perfectly so the max exposure time should
+    # be equal to the individual time times two.
+    time_difference = (
+        output_model.meta.resample.product_exposure_time
+        - 2 * exposure_1[0].meta.exposure.effective_exposure_time
     )
+    assert np.abs(time_difference) < 0.1
 
     assert output_model.meta.exposure.start_time == min(
         x.meta.exposure.start_time for x in input_models
@@ -502,9 +515,13 @@ def test_update_exposure_times_same_sca_different_exposures(exposure_1, exposure
         x.meta.exposure.end_time for x in input_models
     )
 
-    assert output_model.meta.resample.product_exposure_time == sum(
-        x.meta.exposure.exposure_time for x in input_models
+    # likewise the per-pixel median exposure time is just 2x the individual
+    # sca exposure time.
+    time_difference = (
+        output_model.meta.exposure.exposure_time
+        - 2 * exposure_1[0].meta.exposure.effective_exposure_time
     )
+    assert np.abs(time_difference) < 0.1
 
 
 @pytest.mark.parametrize(