From 9f1bbbd67db92ea588453911be3c0b116de56cf9 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Thu, 19 Dec 2024 15:46:46 -0500 Subject: [PATCH] AL-852: GWCS inverse transform should respect its bounding box (#8554) Co-authored-by: Ned Molter Co-authored-by: Ned Molter Co-authored-by: Melanie Clarke Co-authored-by: Tyler Pauly --- changes/8554.assign_wcs.rst | 2 + changes/8554.resample.rst | 1 + changes/8554.skymatch.rst | 1 + jwst/assign_wcs/miri.py | 2 + jwst/assign_wcs/util.py | 2 +- .../tests/test_outlier_detection.py | 50 +++++++++++-------- jwst/resample/resample_utils.py | 9 +++- jwst/resample/tests/test_resample_step.py | 7 ++- jwst/resample/tests/test_utils.py | 4 +- jwst/skymatch/skyimage.py | 2 +- jwst/tweakreg/tests/test_multichip_jwst.py | 6 +-- pyproject.toml | 14 +++--- 12 files changed, 63 insertions(+), 37 deletions(-) create mode 100644 changes/8554.assign_wcs.rst create mode 100644 changes/8554.resample.rst create mode 100644 changes/8554.skymatch.rst diff --git a/changes/8554.assign_wcs.rst b/changes/8554.assign_wcs.rst new file mode 100644 index 0000000000..c36dd9d167 --- /dev/null +++ b/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. diff --git a/changes/8554.resample.rst b/changes/8554.resample.rst new file mode 100644 index 0000000000..bda0397ae8 --- /dev/null +++ b/changes/8554.resample.rst @@ -0,0 +1 @@ +Ignore the bounding_box in the inverse WCS transform in reproject. diff --git a/changes/8554.skymatch.rst b/changes/8554.skymatch.rst new file mode 100644 index 0000000000..1561af36b6 --- /dev/null +++ b/changes/8554.skymatch.rst @@ -0,0 +1 @@ +Ignore the bounding_box of an observation when computing sky statistics. diff --git a/jwst/assign_wcs/miri.py b/jwst/assign_wcs/miri.py index 59ec01d959..f15938ce88 100644 --- a/jwst/assign_wcs/miri.py +++ b/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 diff --git a/jwst/assign_wcs/util.py b/jwst/assign_wcs/util.py index a577e8a8e4..3735a45bee 100644 --- a/jwst/assign_wcs/util.py +++ b/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] diff --git a/jwst/outlier_detection/tests/test_outlier_detection.py b/jwst/outlier_detection/tests/test_outlier_detection.py index 008ae19c64..7e03300422 100644 --- a/jwst/outlier_detection/tests/test_outlier_detection.py +++ b/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,11 +144,29 @@ 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 @@ -160,11 +174,7 @@ def we_many_sci( # 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): diff --git a/jwst/resample/resample_utils.py b/jwst/resample/resample_utils.py index 4d126ed28c..da917d7e41 100644 --- a/jwst/resample/resample_utils.py +++ b/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 = [] diff --git a/jwst/resample/tests/test_resample_step.py b/jwst/resample/tests/test_resample_step.py index 75eddb8530..c9eb4ec1b5 100644 --- a/jwst/resample/tests/test_resample_step.py +++ b/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])) diff --git a/jwst/resample/tests/test_utils.py b/jwst/resample/tests/test_utils.py index 7e7deb3070..518fc61998 100644 --- a/jwst/resample/tests/test_utils.py +++ b/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) diff --git a/jwst/skymatch/skyimage.py b/jwst/skymatch/skyimage.py index 65bc943bec..8a6ca6f3c8 100644 --- a/jwst/skymatch/skyimage.py +++ b/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) diff --git a/jwst/tweakreg/tests/test_multichip_jwst.py b/jwst/tweakreg/tests/test_multichip_jwst.py index 4402f7b897..352d089e5a 100644 --- a/jwst/tweakreg/tests/test_multichip_jwst.py +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index a944d7ba44..2d7907ed2d 100644 --- a/pyproject.toml +++ b/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",