wfss_contam use datamodel in with context

When the DataModel is cleaned up (__del__ is called)
during garbage collection it will close any opened
hdulists. This will make any subsequent attempts to load
data from the hdulist fail (with ValueError: I/O operation on closed file).

DataModel prior to the below PR contained a self reference
which created a reference cycle which made python fail to collect
the object during generation 0 passes with the garbage collector.
This meant that the hdulist was not immediately closed when
the model fell out of scope.

Code in wfss_contam relied on this persistance as the model created
in the following line is no longer referenced on subsequent lines.

dimage = datamodels.open(dir_image_name).data

python will attempt to clean up the opened data model and fail,
keeping the hdulist (stored at dir_image_name) open. This allowed
later lines that reference dimage to load the array data within the
hdulist.

Fixing the reference cycle in stdatamodels in this PR:
spacetelescope/stdatamodels#109
makes the datamodel more easily garbage collected which causes
wfss_contam to fail due to the hdulist closing when the model
is collected.

To keep the model in scope, the above line is replaced with
a with context to make it explicit that the datamodel is held
open which it's contents are read.

update changes

CHANGES.rst

@@ -6,6 +6,11 @@ documentation
 
 - Remove references to pub server [#7421]
 
+wfss_contam
+-----------
+
+- Open image models in with context to keep models open while accessing contents [#7425]]
+
 1.9.1 (2023-01-03)
 ==================
 

jwst/wfss_contam/observations.py

@@ -117,36 +117,37 @@ def create_pixel_list(self):
         for dir_image_name in self.dir_image_names:
 
             log.info(f"Using direct image {dir_image_name}")
-            dimage = datamodels.open(dir_image_name).data
-
-            if self.sed_file is None:
-                # Default pipeline will use sed_file=None, so we need to compute
-                # photometry values that used to come from HST-style header keywords.
-                # Set pivlam, in units of microns, based on filter name.
-                pivlam = float(self.filter[1:4]) / 100.
-
-                # Use pixel fluxes from the direct image.
-                self.fluxes[pivlam] = []
-                for i in range(len(self.IDs)):
-                    # This loads lists of pixel flux values for each source
-                    # from the direct image
-                    self.fluxes[pivlam].append(dimage[self.ys[i], self.xs[i]])
-
-            else:
-                # Use an SED file. Need to normalize the object stamps.
-                for ID in self.IDs:
-                    vg = self.seg == ID
-                    dnew = dimage
-                    if self.renormalize:
-                        sum_seg = np.sum(dimage[vg])  # But normalize by the whole flux
-                        if sum_seg != 0:
-                            dimage[vg] /= sum_seg
-                    else:
-                        log.debug("not renormalizing sources to unity")
-
-                self.fluxes["sed"] = []
-                for i in range(len(self.IDs)):
-                    self.fluxes["sed"].append(dnew[self.ys[i], self.xs[i]])
+            with datamodels.open(dir_image_name) as model:
+                dimage = model.data
+
+                if self.sed_file is None:
+                    # Default pipeline will use sed_file=None, so we need to compute
+                    # photometry values that used to come from HST-style header keywords.
+                    # Set pivlam, in units of microns, based on filter name.
+                    pivlam = float(self.filter[1:4]) / 100.
+
+                    # Use pixel fluxes from the direct image.
+                    self.fluxes[pivlam] = []
+                    for i in range(len(self.IDs)):
+                        # This loads lists of pixel flux values for each source
+                        # from the direct image
+                        self.fluxes[pivlam].append(dimage[self.ys[i], self.xs[i]])
+
+                else:
+                    # Use an SED file. Need to normalize the object stamps.
+                    for ID in self.IDs:
+                        vg = self.seg == ID
+                        dnew = dimage
+                        if self.renormalize:
+                            sum_seg = np.sum(dimage[vg])  # But normalize by the whole flux
+                            if sum_seg != 0:
+                                dimage[vg] /= sum_seg
+                        else:
+                            log.debug("not renormalizing sources to unity")
+
+                    self.fluxes["sed"] = []
+                    for i in range(len(self.IDs)):
+                        self.fluxes["sed"].append(dnew[self.ys[i], self.xs[i]])
 
     def disperse_all(self, order, wmin, wmax, sens_waves, sens_resp, cache=False):
         """