JP-1829: Fix resample_spec output size from input images crossing RA=0

CHANGES.rst

@@ -107,14 +107,18 @@ ramp_fitting
 refpix
 ------
 
-- Added code to handle NIR subarrays that use 4 readout amplifiers.  Uses and applies reference pixel signal from
-  available amplifiers and side reference pixel regions, including odd-even column separation if requested [#5926]
+- Added code to handle NIR subarrays that use 4 readout amplifiers.  Uses and
+  applies reference pixel signal from available amplifiers and side reference
+  pixel regions, including odd-even column separation if requested [#5926]
 
-- Fixed a bug introduced in #5926 that affected refpix calibration of 1-amp NIR subarrays [#5937]
+- Fixed a bug introduced in #5926 that affected refpix calibration of 1-amp NIR
+  subarrays [#5937]
 
 resample
 --------
 
+- Fix ``resample_spec`` output size from input images crossing RA=0 [#5929]
+
 - Propagate variance arrays into ``SlitModel`` used as input for ``ResampleSpecStep`` [#5941]
 
 source_catalog

jwst/assign_mtwcs/tests/test_mtwcs.py

@@ -9,10 +9,10 @@
 
 def test_mt_multislit():
     file_path = os.path.join(data.__path__[0], 'test_mt_asn.json')
-    model = datamodels.open(file_path)
-    assert model[0].slits[0].meta.wcs.output_frame.name == 'world'
-    step = AssignMTWcsStep()
-    result = step.run(model)
+    with datamodels.open(file_path) as model:
+        assert model[0].slits[0].meta.wcs.output_frame.name == 'world'
+        step = AssignMTWcsStep()
+        result = step.run(model)
     assert isinstance(result, datamodels.ModelContainer)
     assert len(result[0].slits) == 1
     assert result[0].slits[0].meta.wcs.output_frame.name == 'moving_target'

jwst/outlier_detection/tests/test_outlier_detection.py

@@ -28,6 +28,7 @@ def sci_blot_image_pair():
 
     sci.data = np.random.normal(loc=background, size=shape, scale=sigma)
     sci.err = np.zeros(shape) + sigma
+    sci.var_rnoise += 0
 
     # Add a source in the center
     sci.data[10, 10] += 20 * sigma
@@ -113,6 +114,7 @@ def we_three_sci():
     sci1.meta.wcs = create_fitswcs(sci1)
 
     sci1.err = np.zeros(shape) + sigma
+    sci1.var_rnoise += 0
 
     # Make copies with different noise
     sci2 = sci1.copy()

jwst/regtest/test_nircam_image.py

@@ -58,8 +58,7 @@ def run_image3pipeline(run_image2pipeline, rtdata_module, jail):
     ]
     for rate_file in rate_files:
         rtdata.get_data(rate_file)
-        args = ["config/calwebb_image2.cfg", rtdata.input,
-                "--steps.resample.skip=True"]
+        args = ["config/calwebb_image2.cfg", rtdata.input]
         Step.from_cmdline(args)
 
     # Get the level3 assocation json file (though not its members) and run

jwst/resample/resample_spec.py

@@ -142,20 +142,20 @@ def build_interpolated_output_wcs(self, refmodel=None):
             # first input model sets intializes wavelength array and defines
             # the spatial scale of the output wcs
             if im == 0:
-                for iw in wavelength_array:
-                    all_wavelength.append(iw)
+                all_wavelength = np.append(all_wavelength, wavelength_array)
 
+                # find the center ra and dec for this slit at central wavelength
                 lam_center_index = int((bb[spectral_axis][1] -
                                         bb[spectral_axis][0]) / 2)
-                if spatial_axis == 0:
-                    ra_center = ra[lam_center_index, :]
-                    dec_center = dec[lam_center_index, :]
+                if spatial_axis == 0:  # MIRI LRS, the WCS x axis is spatial
+                    ra_slice = ra[lam_center_index, :]
+                    dec_slice = dec[lam_center_index, :]
                 else:
-                    ra_center = ra[:, lam_center_index]
-                    dec_center = dec[:, lam_center_index]
-                # find the ra and dec for this slit using center of slit
-                ra_center_pt = np.nanmean(ra_center)
-                dec_center_pt = np.nanmean(dec_center)
+                    ra_slice = ra[:, lam_center_index]
+                    dec_slice = dec[:, lam_center_index]
+                # wrap RA if near zero
+                ra_center_pt = np.nanmean(wrap_ra(ra_slice))
+                dec_center_pt = np.nanmean(dec_slice)
 
                 if resample_utils.is_sky_like(model.meta.wcs.output_frame):
                     # convert ra and dec to tangent projection
@@ -173,7 +173,7 @@ def build_interpolated_output_wcs(self, refmodel=None):
                     x_tan, y_tan = ra, dec
 
                 # pull out data from center
-                if spectral_axis == 0:
+                if spectral_axis == 0:  # MIRI LRS, the WCS x axis is spatial
                     x_tan_array = x_tan.T[lam_center_index]
                     y_tan_array = y_tan.T[lam_center_index]
                 else:
@@ -195,30 +195,24 @@ def build_interpolated_output_wcs(self, refmodel=None):
 
             # append all ra and dec values to use later to find min and max
             # ra and dec
-            ra_use = ra.flatten()
-            ra_use = ra_use[~np.isnan(ra_use)]
-            dec_use = dec.flatten()
-            dec_use = dec_use[~np.isnan(dec_use)]
-            all_ra_slit.append(ra_use)
-            all_dec_slit.append(dec_use)
+            ra_use = ra[~np.isnan(ra)].flatten()
+            dec_use = dec[~np.isnan(dec)].flatten()
+            all_ra_slit = np.append(all_ra_slit, ra_use)
+            all_dec_slit = np.append(all_dec_slit, dec_use)
 
             # now check wavelength array to see if we need to add to it
             this_minw = np.min(wavelength_array)
             this_maxw = np.max(wavelength_array)
             all_minw = np.min(all_wavelength)
             all_maxw = np.max(all_wavelength)
-
             if this_minw < all_minw:
                 addpts = wavelength_array[wavelength_array < all_minw]
-                for ip in range(len(addpts)):
-                    all_wavelength.append(addpts[ip])
+                all_wavelength = np.append(all_wavelength, addpts)
             if this_maxw > all_maxw:
                 addpts = wavelength_array[wavelength_array > all_maxw]
-                for ip in range(len(addpts)):
-                    all_wavelength.append(addpts[ip])
+                all_wavelength = np.append(all_wavelength, addpts)
 
         # done looping over set of models
-
         all_ra = np.hstack(all_ra_slit)
         all_dec = np.hstack(all_dec_slit)
         all_wave = np.hstack(all_wavelength)
@@ -280,11 +274,12 @@ def build_interpolated_output_wcs(self, refmodel=None):
         all_ra = wrap_ra(all_ra)
         ra_min = np.amin(all_ra)
         ra_max = np.amax(all_ra)
-
         ra_center_final = (ra_max + ra_min) / 2.0
+
         dec_min = np.amin(all_dec)
         dec_max = np.amax(all_dec)
         dec_center_final = (dec_max + dec_min) / 2.0
+
         tan = Pix2Sky_TAN()
         if len(self.input_models) == 1:  # single model use ra_center_pt to be consistent
             # with how resample was done before
@@ -311,7 +306,7 @@ def build_interpolated_output_wcs(self, refmodel=None):
         if resample_utils.is_sky_like(model.meta.wcs.output_frame):
             transform = mapping | (pix_to_xtan & pix_to_ytan | undist2sky) & pix_to_wavelength
         else:
-            transform = mapping | (pix_to_xtan & pix_to_ytan) & pix_to_wavelength
+            transform = mapping | pix_to_xtan & pix_to_ytan & pix_to_wavelength
 
         det = cf.Frame2D(name='detector', axes_order=(0, 1))
         if resample_utils.is_sky_like(model.meta.wcs.output_frame):

jwst/resample/resample_spec_step.py

@@ -94,6 +94,7 @@ def _process_multislit(self, input_models):
 
         for container in containers.values():
             resamp = resample_spec.ResampleSpecData(container, **self.drizpars)
+
             drizzled_models = resamp.do_drizzle()
 
             for model in drizzled_models:
@@ -138,6 +139,7 @@ def _process_slit(self, input_models):
 
         resamp = resample_spec.ResampleSpecData(input_models, **self.drizpars)
 
+        # Only drizzle the area within the bounding box
         if input_models[0].meta.exposure.type == "MIR_LRS-FIXEDSLIT":
             bb = input_models[0].meta.wcs.bounding_box
             ((x1, x2), (y1, y2)) = bb

jwst/resample/tests/test_resample_step.py

@@ -7,6 +7,7 @@
 from jwst.assign_wcs import AssignWcsStep
 from jwst.extract_2d import Extract2dStep
 from jwst.resample import ResampleSpecStep, ResampleStep
+from jwst.resample.resample_spec import ResampleSpecData
 
 
 @pytest.fixture
@@ -15,6 +16,7 @@ def nirspec_rate():
     xsize = 2048
     shape = (ysize, xsize)
     im = ImageModel(shape)
+    im.var_rnoise += 1
     im.meta.wcsinfo = {
         'dec_ref': 40,
         'ra_ref': 100,
@@ -71,7 +73,7 @@ def miri_rate():
     shape = (ysize, xsize)
     im = ImageModel(shape)
     im.data += 5
-    im.var_rnoise = np.random.random(shape)
+    im.var_rnoise += 1
     im.meta.wcsinfo = {
         'dec_ref': 40,
         'ra_ref': 100,
@@ -247,6 +249,7 @@ def test_pixel_scale_ratio_imaging(nircam_rate, ratio):
 def test_weight_type(nircam_rate, _jail):
     """Check that weight_type of exptime and ivm work"""
     im1 = AssignWcsStep.call(nircam_rate, sip_approx=False)
+    im1.var_rnoise[:] = 0
     im2 = im1.copy()
     im3 = im1.copy()
     im1.data += 10
@@ -292,3 +295,88 @@ def test_sip_coeffs_do_not_propagate(nircam_rate):
 
     # Make sure we have a PC matrix
     assert result.meta.wcsinfo.pc1_1 is not None
+
+
+@pytest.fixture
+def miri_rate_zero_crossing():
+    xsize = 1032
+    ysize = 1024
+    shape = (ysize, xsize)
+    im = ImageModel(shape)
+    im.var_rnoise = np.random.random(shape)
+    im.meta.wcsinfo = {
+        'dec_ref': 2.16444343946559e-05,
+        'ra_ref': -0.00026031780056776,
+        'roll_ref': 0.0,
+        'v2_ref': -415.0690466121227,
+        'v3_ref': -400.575920398547,
+        'v3yangle': 0.0,
+        'vparity': -1}
+    im.meta.instrument = {
+        'detector': 'MIRIMAGE',
+        'filter': 'P750L',
+        'name': 'MIRI'}
+    im.meta.observation = {
+        'date': '2019-01-01',
+        'time': '17:00:00'}
+    im.meta.subarray = {
+        'fastaxis': 1,
+        'name': 'FULL',
+        'slowaxis': 2,
+        'xsize': xsize,
+        'xstart': 1,
+        'ysize': ysize,
+        'ystart': 1}
+    im.meta.exposure = {
+        'duration': 11.805952,
+        'end_time': 58119.85416,
+        'exposure_time': 11.776,
+        'frame_time': 0.11776,
+        'group_time': 0.11776,
+        'groupgap': 0,
+        'integration_time': 11.776,
+        'nframes': 1,
+        'ngroups': 100,
+        'nints': 1,
+        'nresets_between_ints': 0,
+        'nsamples': 1,
+        'readpatt': 'FAST',
+        'sample_time': 10.0,
+        'start_time': 58119.8333,
+        'type': 'MIR_LRS-FIXEDSLIT',
+        'zero_frame': False}
+
+    return im
+
+
+@pytest.fixture
+def miri_rate_pair(miri_rate_zero_crossing):
+    im1 = miri_rate_zero_crossing
+    # Create a nodded version
+    im2 = im1.copy()
+    im2.meta.wcsinfo.ra_ref = 0.00026308279776455
+    im2.meta.wcsinfo.dec_ref = -2.1860888891293e-05
+    im1 = AssignWcsStep.call(im1)
+    im2 = AssignWcsStep.call(im2)
+
+    return im1, im2
+
+
+def test_build_interpolated_output_wcs(miri_rate_pair):
+    im1, im2 = miri_rate_pair
+
+    driz = ResampleSpecData(ModelContainer([im1, im2]))
+    output_wcs = driz.build_interpolated_output_wcs()
+
+    # Make sure that all RA, Dec values in the input image have a location in
+    # the output frame
+    grid = grid_from_bounding_box(im2.meta.wcs.bounding_box)
+    ra, dec, lam = im2.meta.wcs(*grid)
+    x, y = output_wcs.invert(ra, dec, lam)
+
+    # This currently fails, as we see a slight offset
+    # assert (x > 0).all()
+
+    # Make sure the output slit size is larger than the input slit size
+    # for this nodded data
+    assert driz.data_size[1] > ra.shape[1]