Skip to content

Commit

Permalink
tweakreg into stcal
Browse files Browse the repository at this point in the history
  • Loading branch information
@emolter
emolter committed Jul 3, 2024
1 parent ce9a135 commit fc20b8c
Show file tree
Hide file tree
Showing 12 changed files with 104 additions and 781 deletions.
3 changes: 1 addition & 2 deletions jwst/assign_wcs/__init__.py
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
from .assign_wcs_step import AssignWcsStep
from .nirspec import (nrs_wcs_set_input, nrs_ifu_wcs, get_spectral_order_wrange)
from .niriss import niriss_soss_set_input
from .util import update_fits_wcsinfo

__all__ = ['AssignWcsStep', "nrs_wcs_set_input", "nrs_ifu_wcs", "get_spectral_order_wrange",
"niriss_soss_set_input", "update_fits_wcsinfo"]
"niriss_soss_set_input"]
4 changes: 2 additions & 2 deletions jwst/assign_wcs/assign_wcs_step.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,8 @@
from ..lib.exposure_types import IMAGING_TYPES
import logging
from .assign_wcs import load_wcs
from .util import MSAFileError, update_fits_wcsinfo
from .util import wfss_imaging_wcs, wcs_bbox_from_shape
from .util import MSAFileError
from .util import wfss_imaging_wcs, wcs_bbox_from_shape, update_fits_wcsinfo
from .nircam import imaging as nircam_imaging
from .niriss import imaging as niriss_imaging

Expand Down
164 changes: 2 additions & 162 deletions jwst/assign_wcs/util.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,9 +34,9 @@
_MAX_SIP_DEGREE = 6


__all__ = ["reproject", "wcs_from_footprints", "velocity_correction",
__all__ = ["reproject", "velocity_correction",
"MSAFileError", "NoDataOnDetectorError", "compute_scale",
"calc_rotation_matrix", "wrap_ra", "update_fits_wcsinfo"]
"calc_rotation_matrix", "wrap_ra", "update_fits_wcsinfo",]


class MSAFileError(Exception):
Expand Down Expand Up @@ -197,164 +197,6 @@ def calc_rotation_matrix(roll_ref: float, v3i_yang: float, vparity: int = 1) ->
return [pc1_1, pc1_2, pc2_1, pc2_2]


def wcs_from_footprints(dmodels, refmodel=None, transform=None, bounding_box=None,
pscale_ratio=None, pscale=None, rotation=None,
shape=None, crpix=None, crval=None, wcslist=None):
"""
Create a WCS from a list of input data models.
A fiducial point in the output coordinate frame is created from the
footprints of all WCS objects. For a spatial frame this is the center
of the union of the footprints. For a spectral frame the fiducial is in
the beginning of the footprint range.
If ``refmodel`` is None, the first WCS object in the list is considered
a reference. The output coordinate frame and projection (for celestial frames)
is taken from ``refmodel``.
If ``transform`` is not supplied, a compound transform is created using
CDELTs and PC.
If ``bounding_box`` is not supplied, the bounding_box of the new WCS is computed
from bounding_box of all input WCSs.
Parameters
----------
dmodels : list of `~jwst.datamodels.JwstDataModel`
A list of data models.
refmodel : `~jwst.datamodels.JwstDataModel`, optional
This model's WCS is used as a reference.
WCS. The output coordinate frame, the projection and a
scaling and rotation transform is created from it. If not supplied
the first model in the list is used as ``refmodel``.
transform : `~astropy.modeling.core.Model`, optional
A transform, passed to :meth:`~gwcs.wcstools.wcs_from_fiducial`
If not supplied Scaling | Rotation is computed from ``refmodel``.
bounding_box : tuple, optional
Bounding_box of the new WCS.
If not supplied it is computed from the bounding_box of all inputs.
pscale_ratio : float, None, optional
Ratio of input to output pixel scale. Ignored when either
``transform`` or ``pscale`` are provided.
pscale : float, None, optional
Absolute pixel scale in degrees. When provided, overrides
``pscale_ratio``. Ignored when ``transform`` is provided.
rotation : float, None, optional
Position angle of output image’s Y-axis relative to North.
A value of 0.0 would orient the final output image to be North up.
The default of `None` specifies that the images will not be rotated,
but will instead be resampled in the default orientation for the camera
with the x and y axes of the resampled image corresponding
approximately to the detector axes. Ignored when ``transform`` is
provided.
shape : tuple of int, None, optional
Shape of the image (data array) using ``numpy.ndarray`` convention
(``ny`` first and ``nx`` second). This value will be assigned to
``pixel_shape`` and ``array_shape`` properties of the returned
WCS object.
crpix : tuple of float, None, optional
Position of the reference pixel in the image array. If ``crpix`` is not
specified, it will be set to the center of the bounding box of the
returned WCS object.
crval : tuple of float, None, optional
Right ascension and declination of the reference pixel. Automatically
computed if not provided.
"""
bb = bounding_box
if wcslist is None:
wcslist = [im.meta.wcs for im in dmodels]

if not isiterable(wcslist):
raise ValueError("Expected 'wcslist' to be an iterable of WCS objects.")

if not all([isinstance(w, WCS) for w in wcslist]):
raise TypeError("All items in wcslist are to be instances of gwcs.WCS.")

if refmodel is None:
refmodel = dmodels[0]
else:
if not isinstance(refmodel, JwstDataModel):
raise TypeError("Expected refmodel to be an instance of DataModel.")

fiducial = compute_fiducial(wcslist, bb)
if crval is not None:
# overwrite spatial axes with user-provided CRVAL:
i = 0
for k, axt in enumerate(wcslist[0].output_frame.axes_type):
if axt == 'SPATIAL':
fiducial[k] = crval[i]
i += 1

ref_fiducial = np.array([refmodel.meta.wcsinfo.ra_ref, refmodel.meta.wcsinfo.dec_ref])

prj = astmodels.Pix2Sky_TAN()

if transform is None:
transform = []
wcsinfo = pointing.wcsinfo_from_model(refmodel)
sky_axes, spec, other = gwutils.get_axes(wcsinfo)

# Need to put the rotation matrix (List[float, float, float, float])
# returned from calc_rotation_matrix into the correct shape for
# constructing the transformation
v3yangle = np.deg2rad(refmodel.meta.wcsinfo.v3yangle)
vparity = refmodel.meta.wcsinfo.vparity
if rotation is None:
roll_ref = np.deg2rad(refmodel.meta.wcsinfo.roll_ref)
else:
roll_ref = np.deg2rad(rotation) + (vparity * v3yangle)

pc = np.reshape(
calc_rotation_matrix(roll_ref, v3yangle, vparity=vparity),
(2, 2)
)

rotation = astmodels.AffineTransformation2D(pc, name='pc_rotation_matrix')
transform.append(rotation)

if sky_axes:
if not pscale:
pscale = compute_scale(refmodel.meta.wcs, ref_fiducial,
pscale_ratio=pscale_ratio)
transform.append(astmodels.Scale(pscale, name='cdelt1') & astmodels.Scale(pscale, name='cdelt2'))

if transform:
transform = functools.reduce(lambda x, y: x | y, transform)

out_frame = refmodel.meta.wcs.output_frame
input_frame = refmodel.meta.wcs.input_frame
wnew = wcs_from_fiducial(fiducial, coordinate_frame=out_frame, projection=prj,
transform=transform, input_frame=input_frame)

footprints = [w.footprint().T for w in wcslist]
domain_bounds = np.hstack([wnew.backward_transform(*f) for f in footprints])
axis_min_values = np.min(domain_bounds, axis=1)
domain_bounds = (domain_bounds.T - axis_min_values).T

output_bounding_box = []
for axis in out_frame.axes_order:
axis_min, axis_max = domain_bounds[axis].min(), domain_bounds[axis].max()
output_bounding_box.append((axis_min, axis_max))

output_bounding_box = tuple(output_bounding_box)
if crpix is None:
offset1, offset2 = wnew.backward_transform(*fiducial)
offset1 -= axis_min_values[0]
offset2 -= axis_min_values[1]
else:
offset1, offset2 = crpix
offsets = astmodels.Shift(-offset1, name='crpix1') & astmodels.Shift(-offset2, name='crpix2')

wnew.insert_transform('detector', offsets, after=True)
wnew.bounding_box = output_bounding_box

if shape is None:
shape = [int(axs[1] - axs[0] + 0.5) for axs in output_bounding_box[::-1]]

wnew.pixel_shape = shape[::-1]
wnew.array_shape = shape

return wnew


def compute_fiducial(wcslist, bounding_box=None):
"""
For a celestial footprint this is the center.
Expand Down Expand Up @@ -1283,7 +1125,6 @@ def update_fits_wcsinfo(datamodel, max_pix_error=0.01, degree=None,
Parameters
----------
datamodel : `ImageModel`
The input data model for imaging or WFSS mode whose ``meta.wcsinfo``
field should be updated from GWCS. By default, ``datamodel.meta.wcs``
Expand Down Expand Up @@ -1386,7 +1227,6 @@ def update_fits_wcsinfo(datamodel, max_pix_error=0.01, degree=None,
Notes
-----
Use of this requires a judicious choice of required accuracies.
Attempts to use higher degrees (~7 or higher) will typically fail due
to floating point problems that arise with high powers.
Expand Down
2 changes: 1 addition & 1 deletion jwst/resample/resample.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
from jwst.datamodels import ModelContainer

from . import gwcs_drizzle
from . import resample_utils
from jwst.resample import resample_utils
from ..model_blender import blendmeta

log = logging.getLogger(__name__)
Expand Down
31 changes: 18 additions & 13 deletions jwst/resample/resample_utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,8 @@
from stdatamodels.dqflags import interpret_bit_flags
from stdatamodels.jwst.datamodels.dqflags import pixel

from jwst.assign_wcs.util import wcs_from_footprints, wcs_bbox_from_shape
from jwst.assign_wcs.util import wcs_bbox_from_shape
from stcal.alignment import util


log = logging.getLogger(__name__)
Expand All @@ -22,21 +23,23 @@
def make_output_wcs(input_models, ref_wcs=None,
pscale_ratio=None, pscale=None, rotation=None, shape=None,
crpix=None, crval=None):
""" Generate output WCS here based on footprints of all input WCS objects
"""Generate output WCS here based on footprints of all input WCS objects.
Parameters
----------
input_models : list of `~jwst.datamodel.JwstDataModel`
input_models : list of `DataModel objects`
Each datamodel must have a ~gwcs.WCS object.
pscale_ratio : float, optional
Ratio of input to output pixel scale. Ignored when ``pscale`` is provided.
Ratio of input to output pixel scale. Ignored when ``pscale``
is provided.
pscale : float, None, optional
Absolute pixel scale in degrees. When provided, overrides
``pscale_ratio``.
rotation : float, None, optional
Position angle of output image’s Y-axis relative to North.
Position angle of output image Y-axis relative to North.
A value of 0.0 would orient the final output image to be North up.
The default of `None` specifies that the images will not be rotated,
but will instead be resampled in the default orientation for the camera
Expand All @@ -50,7 +53,7 @@ def make_output_wcs(input_models, ref_wcs=None,
WCS object.
crpix : tuple of float, None, optional
Position of the reference pixel in the image array. If ``crpix`` is not
Position of the reference pixel in the image array. If ``crpix`` is not
specified, it will be set to the center of the bounding box of the
returned WCS object.
Expand All @@ -62,7 +65,6 @@ def make_output_wcs(input_models, ref_wcs=None,
-------
output_wcs : object
WCS object, with defined domain, covering entire set of input frames
"""
if ref_wcs is None:
wcslist = [i.meta.wcs for i in input_models]
Expand All @@ -72,10 +74,11 @@ def make_output_wcs(input_models, ref_wcs=None,
naxes = wcslist[0].output_frame.naxes

if naxes != 2:
raise RuntimeError("Output WCS needs 2 spatial axes. "
f"{wcslist[0]} has {naxes}.")
msg = f"Output WCS needs 2 spatial axes \
but the supplied WCS has {naxes} axes."
raise RuntimeError(msg)

output_wcs = wcs_from_footprints(
output_wcs = util.wcs_from_footprints(
input_models,
pscale_ratio=pscale_ratio,
pscale=pscale,
Expand All @@ -88,15 +91,17 @@ def make_output_wcs(input_models, ref_wcs=None,
else:
naxes = ref_wcs.output_frame.naxes
if naxes != 2:
raise RuntimeError("Output WCS needs 2 spatial axes but the "
f"supplied WCS has {naxes} axes.")
msg = f"Output WCS needs 2 spatial axes \
but the supplied WCS has {naxes} axes."
raise RuntimeError(msg)
output_wcs = deepcopy(ref_wcs)
if shape is not None:
output_wcs.array_shape = shape

# Check that the output data shape has no zero length dimensions
if not np.prod(output_wcs.array_shape):
raise ValueError(f"Invalid output frame shape: {tuple(output_wcs.array_shape)}")
msg = f"Invalid output frame shape: {tuple(output_wcs.array_shape)}"
raise ValueError(msg)

return output_wcs

Expand Down
Loading

0 comments on commit fc20b8c

Please sign in to comment.