Repository: spacetelescope/jwst

Issue:
Outlier blotting does not account for img.meta.photometry.pixelarea_steradians
_Issue [JP-3635](https://jira.stsci.edu/browse/JP-3635) was created on JIRA by [David Law](https://jira.stsci.edu/secure/ViewProfile.jspa?name=dlaw):_

As discussed by <https://jira.stsci.edu/browse/JP-3548,> a small error in the MIRI PHOTOM/AREA files exposed an error in the outlier detection blotting routine.

In brief, resample accounts for the difference between the nominal pixel area in steradians used for the photometric calibration (img.meta.photometry.pixelarea_steradians) and the as-computed pixel area from the WCS when resampling cal files.  That isn't performed in reverse though when blotting median images back to match those individual cal files.

See <https://stsci.box.com/s/8nn6nqo0i244wzhtnelezqmabdargtv1> for an example test case.  The cal files all get about 2% fainter when mapped to i2d and median files, and stay faint when mapped back to blot files.  As a result, it's very easy to exceed the matching tolerance and flag most of the image as cosmic rays.

Discussion:
emolter: Regression tests [started here](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1518/).
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/spacetelescope/jwst/pull/8553?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 58.54%. Comparing base [(`b7e0b10`)](https://app.codecov.io/gh/spacetelescope/jwst/commit/b7e0b107d54ee411b252518e149bbf5501e10403?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) to head [(`1ba3852`)](https://app.codecov.io/gh/spacetelescope/jwst/commit/1ba3852fdc44091be6a97a59228b2cd5771f8e58?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).
> Report is 9 commits behind head on master.


<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #8553      +/-   ##
==========================================
+ Coverage   58.02%   58.54%   +0.52%     
==========================================
  Files         388      388              
  Lines       38977    38960      -17     
==========================================
+ Hits        22617    22810     +193     
+ Misses      16360    16150     -210     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/spacetelescope/jwst/pull/8553?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).

braingram: @mairanteodoro Does romancal need a similar fix?
drlaw1558: Looks good to me; as a test I tried deliberately messing up the header PIXAR keywords for some MIRI data by 20% to really throw off the outlier detection routine.  In the image attached the left panel is the original drizzled i2d image (with good PIXAR values), the middle panel is the existing jwst code working on the data with bad PIXAR values, and the right panel is this PR working on the data with bad PIXAR values.  As expected the current code does really bad things because the scaling isn't being applied when blotting.  This PR looks like it cleans that up really well, and the result looks identical to the original (except, of course, for the fluxcal error that I introduced by hacking the PIXAR keywords).

<img width="1558" alt="jp3635" src="https://github.com/spacetelescope/jwst/assets/11397442/ce498bb2-d8eb-48b1-be41-2556b06e81d9">

emolter: Most of the 80 or so regtest failures stem from `ValueError: WCS must have array_shape attribute set.` when attempting to do `compute_image_pixel_area(blot_img.meta.wcs)` in this step.  The modes affected are seemingly all NIRSpec modes as well as MIRI LRS.  Do we expect the wcs to have an `array_shape` for these modes?  

The `array_shape` seems to be set by `resample_step.get_drizpars()` with errors if it cannot be computed.  We only attempt to do blotting if `resample_data is True`.  However, the `outlier_detection` step does not run the entire resample step, calling `resample.ResampleData` directly instead.  So the options seem to be to run `get_drizpars` from inside this step (not sure what complications might arise), to hard-code `pix_ratio` to unity if `array_shape` is not set, or to manually set the `array_shape` (not sure precisely how to do this yet).
melanieclarke: > Most of the 80 or so regtest failures stem from `ValueError: WCS must have array_shape attribute set.` when attempting to do `compute_image_pixel_area(blot_img.meta.wcs)` in this step. The modes affected are seemingly all NIRSpec modes as well as MIRI LRS. Do we expect the wcs to have an `array_shape` for these modes?

Currently, spectral data does not have iscale computed/applied in resampling, so it shouldn't be applied in blotting either.  You'll need to check for `'SPECTRAL' not in img.meta.wcs.output_frame.axes_type`, as is done in the resample module.
emolter: started [another regtest run](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1523/) with the suggested change
hbushouse: @emolter Your latest regtest run is still corrupted with unrelated changes due to 1) new NIRCam photom ref files and 2) changes to the jump step in stcal. The truth files for the NRC photom related tests have been updated now, and a new version of stcal (1.7.2) has been released that contains the jump updates (and truth files updated). So try running again, being sure you're using stcal 1.7.2 (or stcal/main).
emolter: > @emolter Your latest regtest run is still corrupted with unrelated changes due to 1) new NIRCam photom ref files and 2) changes to the jump step in stcal. The truth files for the NRC photom related tests have been updated now, and a new version of stcal (1.7.2) has been released that contains the jump updates (and truth files updated). So try running again, being sure you're using stcal 1.7.2 (or stcal/main).


[started here](https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1526/)



Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 375398f4a..c9ba6302f 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -171,6 +171,9 @@ outlier_detection
   to detect outliers in TSO data, with user-defined
   rolling window width via the ``rolling_window_width`` parameter. [#8473]
 
+- Fixed a bug that led to small total flux offsets between input and blotted
+  images if the nominal and actual wcs-computed pixel areas were different. [#8553]
+
 pathloss
 --------
 
diff --git a/jwst/outlier_detection/outlier_detection.py b/jwst/outlier_detection/outlier_detection.py
index a44a89433..1e5d3c50c 100644
--- a/jwst/outlier_detection/outlier_detection.py
+++ b/jwst/outlier_detection/outlier_detection.py
@@ -17,6 +17,7 @@ from jwst.datamodels import ModelContainer
 from jwst.resample import resample
 from jwst.resample.resample_utils import build_driz_weight, calc_gwcs_pixmap
 from jwst.stpipe import Step
+from jwst.resample.resample import compute_image_pixel_area
 
 log = logging.getLogger(__name__)
 log.setLevel(logging.DEBUG)
@@ -568,7 +569,12 @@ def gwcs_blot(median_model, blot_img, interp='poly5', sinscl=1.0):
     log.debug("Pixmap shape: {}".format(pixmap[:, :, 0].shape))
     log.debug("Sci shape: {}".format(blot_img.data.shape))
 
-    pix_ratio = 1
+    if 'SPECTRAL' not in blot_img.meta.wcs.output_frame.axes_type:
+        input_pixflux_area = blot_img.meta.photometry.pixelarea_steradians
+        input_pixel_area = compute_image_pixel_area(blot_img.meta.wcs)
+        pix_ratio = np.sqrt(input_pixflux_area / input_pixel_area)
+    else:
+        pix_ratio = 1.0
     log.info('Blotting {} <-- {}'.format(blot_img.data.shape, median_model.data.shape))
 
     outsci = np.zeros(blot_img.shape, dtype=np.float32)
diff --git a/jwst/resample/resample.py b/jwst/resample/resample.py
index 725ee94bc..1a9ff734c 100644
--- a/jwst/resample/resample.py
+++ b/jwst/resample/resample.py
@@ -114,7 +114,7 @@ class ResampleData:
                 self.output_wcs.array_shape = output_shape[::-1]
 
             if output_wcs.pixel_area is None:
-                output_pix_area = _compute_image_pixel_area(self.output_wcs)
+                output_pix_area = compute_image_pixel_area(self.output_wcs)
                 if output_pix_area is None:
                     raise ValueError(
                         "Unable to compute output pixel area from 'output_wcs'."
@@ -250,7 +250,7 @@ class ResampleData:
                 if (input_pixflux_area and
                         'SPECTRAL' not in img.meta.wcs.output_frame.axes_type):
                     img.meta.wcs.array_shape = img.data.shape
-                    input_pixel_area = _compute_image_pixel_area(img.meta.wcs)
+                    input_pixel_area = compute_image_pixel_area(img.meta.wcs)
                     if input_pixel_area is None:
                         raise ValueError(
                             "Unable to compute input pixel area from WCS of input "
@@ -336,7 +336,7 @@ class ResampleData:
             if (input_pixflux_area and
                     'SPECTRAL' not in img.meta.wcs.output_frame.axes_type):
                 img.meta.wcs.array_shape = img.data.shape
-                input_pixel_area = _compute_image_pixel_area(img.meta.wcs)
+                input_pixel_area = compute_image_pixel_area(img.meta.wcs)
                 if input_pixel_area is None:
                     raise ValueError(
                         "Unable to compute input pixel area from WCS of input "
@@ -805,7 +805,7 @@ def _get_boundary_points(xmin, xmax, ymin, ymax, dx=None, dy=None, shrink=0):
     return x, y, area, center, b, r, t, l
 
 
-def _compute_image_pixel_area(wcs):
+def compute_image_pixel_area(wcs):
     """ Computes pixel area in steradians.
     """
     if wcs.array_shape is None:


Test patch:
diff --git a/jwst/outlier_detection/tests/test_outlier_detection.py b/jwst/outlier_detection/tests/test_outlier_detection.py
index 7757ce61d..40982cd4a 100644
--- a/jwst/outlier_detection/tests/test_outlier_detection.py
+++ b/jwst/outlier_detection/tests/test_outlier_detection.py
@@ -141,6 +141,10 @@ def we_many_sci(
     sci1.var_rnoise = np.zeros(shape) + 1.0
     sci1.meta.filename = "foo1_cal.fits"
 
+    # add pixel areas
+    sci1.meta.photometry.pixelarea_steradians = 1.0
+    sci1.meta.photometry.pixelarea_arcsecsq = 1.0
+
     # Make copies with different noise
     all_sci = [sci1]
     for i in range(numsci - 1):
diff --git a/jwst/resample/tests/test_resample_step.py b/jwst/resample/tests/test_resample_step.py
index fc9cd101b..eaf6f1540 100644
--- a/jwst/resample/tests/test_resample_step.py
+++ b/jwst/resample/tests/test_resample_step.py
@@ -12,7 +12,7 @@ from jwst.assign_wcs import AssignWcsStep
 from jwst.assign_wcs.util import compute_fiducial, compute_scale
 from jwst.extract_2d import Extract2dStep
 from jwst.resample import ResampleSpecStep, ResampleStep
-from jwst.resample.resample import _compute_image_pixel_area
+from jwst.resample.resample import compute_image_pixel_area
 from jwst.resample.resample_spec import ResampleSpecData
 
 
@@ -28,7 +28,7 @@ def _set_photom_kwd(im):
         bb = ((xmin - 0.5, xmax - 0.5), (ymin - 0.5, ymax - 0.5))
         im.meta.wcs.bounding_box = bb
 
-    mean_pixel_area = _compute_image_pixel_area(im.meta.wcs)
+    mean_pixel_area = compute_image_pixel_area(im.meta.wcs)
     if mean_pixel_area:
         im.meta.photometry.pixelarea_steradians = mean_pixel_area
         im.meta.photometry.pixelarea_arcsecsq = (