JP-3547: Fix MIRI array sizing for multiple inputs and pixel ratio !=…

… 1 (spacetelescope#8727)

CHANGES.rst

@@ -78,7 +78,7 @@ resample_spec
 
 - Updated handling for the ``pixel_scale_ratio`` parameter to apply only to the
   spatial dimension, to match the sense of the parameter application to the
-  documented intent, and to conserve spectral flux when applied. [#8596]
+  documented intent, and to conserve spectral flux when applied. [#8596, #8727]
 
 - Implemented handling for the ``pixel_scale`` parameter, which was previously
   ignored for spectral resampling. [#8596]

jwst/resample/resample_spec.py

@@ -509,6 +509,7 @@ def build_interpolated_output_wcs(self, input_models):
         all_wavelength = []
         all_ra_slit = []
         all_dec_slit = []
+        xstop = 0
 
         for im, model in enumerate(input_models):
             wcs = model.meta.wcs
@@ -700,10 +701,10 @@ def build_interpolated_output_wcs(self, input_models):
         else:
             pix_to_xtan.intercept = -0.5 * (x_size - 1) * pix_to_xtan.slope
 
-        # single model use size of x_tan_array
+        # single model: use size of x_tan_array
         # to be consistent with method before
         if len(input_models) == 1:
-            x_size = len(x_tan_array)
+            x_size = int(np.ceil(xstop))
 
         # define the output wcs
         transform = mapping | (pix_to_xtan & pix_to_ytan | undist2sky) & pix_to_wavelength
@@ -723,7 +724,7 @@ def build_interpolated_output_wcs(self, input_models):
         # compute the output array size in WCS axes order, i.e. (x, y)
         output_array_size = [0, 0]
         output_array_size[spectral_axis] = int(np.ceil(len(wavelength_array)))
-        output_array_size[spatial_axis] = int(np.ceil(x_size * self.pscale_ratio))
+        output_array_size[spatial_axis] = x_size
 
         # turn the size into a numpy shape in (y, x) order
         output_wcs.array_shape = output_array_size[::-1]

jwst/resample/tests/test_resample_step.py

@@ -490,6 +490,72 @@ def test_pixel_scale_ratio_spec_miri(miri_cal, ratio, units):
     result3.close()
 
 
+@pytest.mark.parametrize("units", ["MJy", "MJy/sr"])
+@pytest.mark.parametrize("ratio", [0.7, 1.0, 1.3])
+def test_pixel_scale_ratio_spec_miri_pair(miri_rate_pair, ratio, units):
+    im1, im2 = miri_rate_pair
+    _set_photom_kwd(im1)
+    _set_photom_kwd(im2)
+    im1.meta.bunit_data = units
+    im2.meta.bunit_data = units
+    im1.meta.filename = 'file1.fits'
+    im2.meta.filename = 'file2.fits'
+    im1.data += 1.0
+    im2.data += 1.0
+
+    # Make an input pixel scale equivalent to the specified ratio
+    input_scale = compute_spectral_pixel_scale(im1.meta.wcs, disp_axis=2)
+    pscale = 3600.0 * input_scale / ratio
+
+    result1 = ResampleSpecStep.call([im1, im2])
+    result2 = ResampleSpecStep.call([im1, im2], pixel_scale_ratio=ratio)
+    result3 = ResampleSpecStep.call([im1, im2], pixel_scale=pscale)
+
+    # pixel_scale and pixel_scale_ratio should be equivalent
+    nn = np.isnan(result2.data) | np.isnan(result3.data)
+    assert np.allclose(result2.data[~nn], result3.data[~nn])
+
+    # Check result2 for expected results
+
+    # wavelength size does not change
+    assert result1.data.shape[0] == result2.data.shape[0]
+
+    # spatial dimension is scaled
+    assert np.isclose(result1.data.shape[1], result2.data.shape[1] / ratio, atol=1)
+
+    # data is non-trivial
+    assert np.nansum(result1.data) > 0.0
+    assert np.nansum(result2.data) > 0.0
+
+    # flux is conserved
+    if 'sr' not in units:
+        # flux density conservation: sum over pixels in each row
+        # needs to be about the same, other than the edges
+        # Check the maximum sums, to avoid edges.
+        assert np.allclose(np.max(np.nansum(result1.data, axis=1)),
+                           np.max(np.nansum(result1.data, axis=1)), rtol=0.05)
+    else:
+        # surface brightness conservation: mean values are the same
+        assert np.allclose(np.nanmean(result1.data, axis=1),
+                           np.nanmean(result2.data, axis=1), rtol=0.05,
+                           equal_nan=True)
+
+    # output area is updated either way
+    area1 = result1.meta.photometry.pixelarea_steradians
+    area2 = result2.meta.photometry.pixelarea_steradians
+    area3 = result2.meta.photometry.pixelarea_steradians
+    assert np.isclose(area1 / area2, ratio)
+    assert np.isclose(area1 / area3, ratio)
+
+    assert result1.meta.resample.pixel_scale_ratio == 1.0
+    assert result2.meta.resample.pixel_scale_ratio == ratio
+    assert np.isclose(result3.meta.resample.pixel_scale_ratio, ratio)
+
+    result1.close()
+    result2.close()
+    result3.close()
+
+
 @pytest.mark.parametrize("units", ["MJy", "MJy/sr"])
 @pytest.mark.parametrize("ratio", [0.7, 1.0, 1.3])
 def test_pixel_scale_ratio_spec_nirspec(nirspec_cal, ratio, units):