Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 25 commits into from
Dec 19, 2024
Merged
Show file tree
Hide file tree
Changes from 24 commits
Commits
Show all changes
25 commits
Select commit Hold shift + click to select a range
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions changes/8554.assign_wcs.rst
Original file line number Diff line number Diff line change
@@ -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.
1 change: 1 addition & 0 deletions changes/8554.resample.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
Ignore the bounding_box in the inverse WCS transform in reproject.
1 change: 1 addition & 0 deletions changes/8554.skymatch.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
Ignore the bounding_box of an observation when computing sky statistics.
2 changes: 2 additions & 0 deletions jwst/assign_wcs/miri.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion jwst/assign_wcs/util.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
50 changes: 30 additions & 20 deletions jwst/outlier_detection/tests/test_outlier_detection.py
Original file line number Diff line number Diff line change
Expand Up @@ -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 (
Expand Down Expand Up @@ -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"
Expand All @@ -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
Expand All @@ -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]
Expand Down Expand Up @@ -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):
Expand Down
9 changes: 8 additions & 1 deletion jwst/resample/resample_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
nden marked this conversation as resolved.
Show resolved Hide resolved
# 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)
Copy link
Member

@mcara mcara Dec 18, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought deepcopying a WCS is cpu-intense. Wouldn't it be better to you try-finally to turn it off/on? Still, I think @emolter 's comment about turning this off may be hiding some other issue is relevant. That is, if we do not disable bbox on the inverse, will the custom WCS test (as modified here) still work?

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 = []
Expand Down
7 changes: 5 additions & 2 deletions jwst/resample/tests/test_resample_step.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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]))
Expand Down
4 changes: 2 additions & 2 deletions jwst/resample/tests/test_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion jwst/skymatch/skyimage.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
6 changes: 3 additions & 3 deletions jwst/tweakreg/tests/test_multichip_jwst.py
Original file line number Diff line number Diff line change
Expand Up @@ -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'

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
14 changes: 7 additions & 7 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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 @ git+https://github.com/spacetelescope/gwcs.git@master",
"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",
Expand Down
Loading