AL-852: GWCS inverse transform should respect its bounding box (#8554)

Co-authored-by: Ned Molter <[email protected]>
Co-authored-by: Ned Molter <[email protected]>
Co-authored-by: Melanie Clarke <[email protected]>
Co-authored-by: Tyler Pauly <[email protected]>

changes/8554.assign_wcs.rst

@@ -0,0 +1,2 @@
+Use the range of points in the TabularModel to adjust the bounding_box in a MIRI LRS FIXEDSLIT observation.
+Ignore the bounding_box in running the inverse WCS transform when computing crpix.

changes/8554.resample.rst

@@ -0,0 +1 @@
+Ignore the bounding_box in the inverse WCS transform in reproject.

changes/8554.skymatch.rst

@@ -0,0 +1 @@
+Ignore the bounding_box of an observation when computing sky statistics.

jwst/assign_wcs/miri.py

@@ -310,6 +310,8 @@ def lrs_xytoabl(input_model, reference_files):
     dxmodel = models.Tabular1D(lookup_table=xshiftref, points=ycen_subarray, name='xshiftref',
                                  bounds_error=False, fill_value=np.nan)
 
+    if input_model.meta.exposure.type.lower() == 'mir_lrs-fixedslit':
+        bb_sub = (bb_sub[0], (dxmodel.points[0].min(), dxmodel.points[0].max()))
     # Fit for the wavelength as a function of Y
     # Reverse the vectors so that yinv is increasing (needed for spline fitting function)
     # Spline fit with enforced smoothness

jwst/assign_wcs/util.py

@@ -127,7 +127,7 @@ def compute_scale(wcs: WCS, fiducial: Union[tuple, np.ndarray],
     if spectral and disp_axis is None:
         raise ValueError('If input WCS is spectral, a disp_axis must be given')
 
-    crpix = np.array(wcs.invert(*fiducial))
+    crpix = np.array(wcs.invert(*fiducial, with_bounding_box=False))
 
     delta = np.zeros_like(crpix)
     spatial_idx = np.where(np.array(wcs.output_frame.axes_type) == 'SPATIAL')[0]

jwst/outlier_detection/tests/test_outlier_detection.py

@@ -6,11 +6,9 @@
 
 from gwcs.wcs import WCS
 from stdatamodels.jwst import datamodels
-from stcal.alignment.util import compute_s_region_imaging
 
 from jwst.datamodels import ModelContainer, ModelLibrary
 from jwst.assign_wcs import AssignWcsStep
-from jwst.assign_wcs.pointing import create_fitswcs
 from jwst.outlier_detection import OutlierDetectionStep
 from jwst.outlier_detection.utils import _flag_resampled_model_crs
 from jwst.outlier_detection.outlier_detection_step import (
@@ -108,21 +106,17 @@ def test_flag_cr(sci_blot_image_pair):
     assert sci.dq[10, 10] == datamodels.dqflags.pixel["GOOD"]
 
 
-# not a fixture - now has options
-def we_many_sci(
-    numsci=3, sigma=0.02, background=1.5, signal=7, exptype="MIR_IMAGE", tsovisit=False
-):
-    """Provide numsci science images with different noise but identical source
-    and same background level"""
-    shape = (21, 20)
+def make_sci1(shape):
+    """Needs to be a fixture because we want to change exposure.type
+    in the subsequent tests without rerunning AssignWCS"""
 
     sci1 = datamodels.ImageModel(shape)
 
     # Populate keywords
     sci1.meta.instrument.name = "MIRI"
     sci1.meta.instrument.detector = "MIRIMAGE"
-    sci1.meta.exposure.type = exptype
-    sci1.meta.visit.tsovisit = tsovisit
+    sci1.meta.exposure.type = "MIR_IMAGE"
+    sci1.meta.visit.tsovisit = False
     sci1.meta.observation.date = "2020-01-01"
     sci1.meta.observation.time = "00:00:00"
     sci1.meta.telescope = "JWST"
@@ -135,6 +129,8 @@ def we_many_sci(
     sci1.meta.wcsinfo.roll_ref = 0
     sci1.meta.wcsinfo.ra_ref = 1.5e-5
     sci1.meta.wcsinfo.dec_ref = 1.5e-5
+    sci1.meta.wcsinfo.v2_ref = 0
+    sci1.meta.wcsinfo.v3_ref = 0
     sci1.meta.wcsinfo.v3yangle = 0
     sci1.meta.wcsinfo.vparity = -1
     sci1.meta.wcsinfo.pc1_1 = 1
@@ -148,23 +144,37 @@ def we_many_sci(
     sci1.meta.wcsinfo.cunit1 = "deg"
     sci1.meta.wcsinfo.cunit2 = "deg"
     sci1.meta.background.subtracted = False
-    sci1.meta.background.level = background
+    sci1.meta.background.level = 1.5
+
+    sci1 = AssignWcsStep.call(sci1)
+
+    sci1.meta.filename = "foo1_cal.fits"
+
+    # add pixel areas
+    sci1.meta.photometry.pixelarea_steradians = 1.0
+    sci1.meta.photometry.pixelarea_arcsecsq = 1.0
+    return sci1
+
+
+# not a fixture - now has options
+def we_many_sci(
+    numsci=3, sigma=0.02, background=1.5, signal=7, exptype="MIR_IMAGE", tsovisit=False
+):
+    """Provide numsci science images with different noise but identical source
+    and same background level"""
+    shape = (21,20)
+    sci1 = make_sci1(shape)
+    sci1.meta.exposure.type = exptype
+    sci1.meta.visit.tsovisit = tsovisit
 
-    # Replace the FITS-type WCS with an Identity WCS
-    sci1.meta.wcs = create_fitswcs(sci1)
-    sci1.meta.wcsinfo.s_region = compute_s_region_imaging(sci1.meta.wcs, shape=shape, center=False)
     rng = np.random.default_rng(720)
     sci1.data = rng.normal(loc=background, size=shape, scale=sigma)
     sci1.err = np.zeros(shape) + sigma
     sci1.data[7, 7] += signal
     # update the noise for this source to include the photon/measurement noise
     sci1.err[7, 7] = np.sqrt(sigma ** 2 + signal)
     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]
@@ -650,7 +660,7 @@ def test_drizzle_and_median_with_resample(three_sci_as_asn, tmp_cwd):
         0.7)
 
     assert isinstance(wcs, WCS)
-    assert median.shape == (21,20)
+    assert median.shape == (34,34)
 
     resamp.single = False
     with pytest.raises(ValueError):

jwst/resample/resample_utils.py

@@ -163,11 +163,18 @@ def reproject(wcs1, wcs2):
     """
 
     try:
+        # Here we want to use the WCS API functions so that a Sliced WCS
+        # will work as well. However, the API functions do not accept
+        # keyword arguments and `with_bounding_box=False` cannot be passsed.
+        # We delete the bounding box on a copy of the WCS - yes, inefficient.
         forward_transform = wcs1.pixel_to_world_values
-        backward_transform = wcs2.world_to_pixel_values
+        wcs_no_bbox = deepcopy(wcs2)
+        wcs_no_bbox.bounding_box = None
+        backward_transform = wcs_no_bbox.world_to_pixel_values
     except AttributeError as err:
         raise TypeError("Input should be a WCS") from err
 
+
     def _reproject(x, y):
         sky = forward_transform(x, y)
         flat_sky = []

jwst/resample/tests/test_resample_step.py

@@ -831,8 +831,8 @@ def test_resample_undefined_variance(nircam_rate, shape):
 
 @pytest.mark.parametrize('ratio', [0.7, 1.2])
 @pytest.mark.parametrize('rotation', [0, 15, 135])
-@pytest.mark.parametrize('crpix', [(256, 488), (700, 124)])
-@pytest.mark.parametrize('crval', [(50, 77), (20, -30)])
+@pytest.mark.parametrize('crpix', [(600, 550), (601, 551)])
+@pytest.mark.parametrize('crval', [(22.04, 11.98), (22.041, 11.981)])
 @pytest.mark.parametrize('shape', [(1205, 1100)])
 def test_custom_wcs_resample_imaging(nircam_rate, ratio, rotation, crpix, crval, shape):
     im = AssignWcsStep.call(nircam_rate, sip_approx=False)
@@ -848,6 +848,9 @@ def test_custom_wcs_resample_imaging(nircam_rate, ratio, rotation, crpix, crval,
 
     t = result.meta.wcs.forward_transform
 
+    # make sure results are nontrivial
+    assert not np.all(np.isnan(result.data))
+
     # test rotation
     pc = t['pc_rotation_matrix'].matrix.value
     orientation = np.rad2deg(np.arctan2(pc[0, 1], pc[1, 1]))

jwst/resample/tests/test_utils.py

@@ -38,7 +38,7 @@ def wcs_gwcs():
     crpix = (500.0, 500.0)
     shape = (1000, 1000)
     pscale = 0.06 / 3600
-    
+
     prj = astmodels.Pix2Sky_TAN()
     fiducial = np.array(crval)
 
@@ -193,7 +193,7 @@ def test_reproject(wcs1, wcs2, offset, request):
     wcs1 = request.getfixturevalue(wcs1)
     wcs2 = request.getfixturevalue(wcs2)
     x = np.arange(150, 200)
-    
+
     f = reproject(wcs1, wcs2)
     res = f(x, x)
     assert_allclose(x, res[0] + offset)

jwst/skymatch/skyimage.py

@@ -619,7 +619,7 @@ def calc_sky(self, overlap=None, delta=True):
                     continue
 
                 # set pixels in 'fill_mask' that are inside a polygon to True:
-                x, y = self.wcs_inv(ra, dec)
+                x, y = self.wcs_inv(ra, dec, with_bounding_box=False)
                 poly_vert = list(zip(*[x, y]))
 
                 polygon = region.Polygon(True, poly_vert)

jwst/tweakreg/tests/test_multichip_jwst.py

@@ -81,7 +81,7 @@ def _make_gwcs_wcs(fits_hdr):
              Mapping((1, 2), name='xtyt'))
     c2tan.name = 'Cartesian 3D to TAN'
 
-    tan2c = (Mapping((0, 0, 1), n_inputs=2, name='xtyt2xyz') |
+    tan2c = (Mapping((0, 0, 1), name='xtyt2xyz') |
              (Const1D(1, name='one') & Identity(2, name='I(2D)')))
     tan2c.name = 'TAN to cartesian 3D'
 
@@ -376,7 +376,7 @@ def test_multichip_alignment_step_rel(monkeypatch):
         format='ascii.ecsv', delimiter=' ',
         names=['RA', 'DEC']
     )
-    x, y = wr.world_to_pixel(refcat['RA'], refcat['DEC'])
+    x, y = wr.invert(refcat['RA'].value, refcat['DEC'].value,  with_bounding_box=False)
     refcat['x'] = x
     refcat['y'] = y
     mr.tweakreg_catalog = refcat
@@ -459,7 +459,7 @@ def test_multichip_alignment_step_abs(monkeypatch):
         format='ascii.ecsv', delimiter=' ',
         names=['RA', 'DEC']
     )
-    x, y = wr.world_to_pixel(refcat['RA'], refcat['DEC'])
+    x, y = wr.invert(refcat['RA'].value, refcat['DEC'].value, with_bounding_box=False)
     refcat['x'] = x
     refcat['y'] = y
     mr.tweakreg_catalog = refcat

pyproject.toml

@@ -18,28 +18,28 @@ classifiers = [
     "Programming Language :: Python :: 3.12",
 ]
 dependencies = [
-    "asdf>=3.1.0,<5",
-    "astropy>=5.3",
+    "asdf>=3.3,<5",
+    "astropy>=6.1",
     "BayesicFitting>=3.0.1",
     "crds>=12.0.3",
     "drizzle>=2.0.0",
-    "gwcs>=0.21.0,<0.23.0",
-    "numpy>=1.22,<2.0",
+    "gwcs>=0.22.0,<0.23.0",
+    "numpy>=1.24,<2.0",
     "opencv-python-headless>=4.6.0.66",
     "photutils>=1.5.0",
     "poppy>=1.0.2",
     "pyparsing>=2.2.1",
     "requests>=2.22",
     "scikit-image>=0.19",
-    "scipy>=1.9.3",
+    "scipy>=1.14.1",
     "spherical-geometry>=1.2.22",
-    "stcal @ git+https://github.com/spacetelescope/stcal.git@30dd76949b6c46ed1ea2a410864252ce6f120c89",
+    "stcal @ git+https://github.com/spacetelescope/stcal.git@main",
     "stdatamodels @ git+https://github.com/spacetelescope/stdatamodels.git@main",
     "stpipe @ git+https://github.com/spacetelescope/stpipe.git@main",
     "stsci.imagestats>=1.6.3",
     "synphot>=1.2",
     "tweakwcs>=0.8.8",
-    "asdf-astropy>=0.3.0",
+    "asdf-astropy>=0.5.0",
     "wiimatch>=0.2.1",
     "packaging>20.0",
     "importlib-metadata>=4.11.4",