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

Switch resample to the new drizzle package API #8866

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
1 change: 1 addition & 0 deletions changes/8866.resample.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
Updated resample code to support the new ``drizzle`` API, see https://github.com/spacetelescope/drizzle/pull/134 for more details.
4 changes: 1 addition & 3 deletions jwst/resample/gwcs_drizzle.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
import numpy as np

from drizzle import util
from drizzle import cdrizzle
from . import resample_utils

Expand Down Expand Up @@ -338,7 +337,7 @@ def dodrizzle(insci, input_wcs, inwht, output_wcs, outsci, outwht, outcon,
"""

# Insure that the fillval parameter gets properly interpreted for use with tdriz
if util.is_blank(str(fillval)):
if str(fillval).strip() == '':
fillval = 'NAN'
else:
fillval = str(fillval)
Expand Down Expand Up @@ -382,7 +381,6 @@ def dodrizzle(insci, input_wcs, inwht, output_wcs, outsci, outwht, outcon,
# Compute the mapping between the input and output pixel coordinates
# for use in drizzle.cdrizzle.tdriz
pixmap = resample_utils.calc_gwcs_pixmap(input_wcs, output_wcs, insci.shape)
# inwht[np.isnan(pixmap[:,:,0])] = 0.

log.debug(f"Pixmap shape: {pixmap[:,:,0].shape}")
log.debug(f"Input Sci shape: {insci.shape}")
Expand Down
49 changes: 27 additions & 22 deletions jwst/resample/resample.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@

import numpy as np
import psutil
from drizzle import cdrizzle, util
from drizzle import cdrizzle
from spherical_geometry.polygon import SphericalPolygon

from stdatamodels.jwst import datamodels
Expand Down Expand Up @@ -192,7 +192,6 @@
output_pix_area * np.rad2deg(3600)**2
)


def do_drizzle(self, input_models):
"""Pick the correct drizzling mode based on self.single
"""
Expand Down Expand Up @@ -251,10 +250,10 @@
else:
iscale = 1.0
return iscale

def resample_group(self, input_models, indices):
"""Apply resample_many_to_many for one group

Parameters
----------
input_models : ModelLibrary
Expand Down Expand Up @@ -290,7 +289,7 @@
if isinstance(img, datamodels.SlitModel):
# must call this explicitly to populate area extension
# although the existence of this extension may not be necessary
img.area = img.area
img.area = img.area

Check warning on line 292 in jwst/resample/resample.py

View check run for this annotation

Codecov / codecov/patch

jwst/resample/resample.py#L292

Added line #L292 was not covered by tests
Copy link
Collaborator

Choose a reason for hiding this comment

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

what does this line do? does it break anything to remove it?

Copy link
Member Author

Choose a reason for hiding this comment

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

I have no idea. It is a strange line but in Python, sometimes things are magical. Seems to be related to SlitModel and also used (set) in the photom step (multi-slit). Maybe @melanieclarke knows more.

iscale = self._get_intensity_scale(img)
log.debug(f'Using intensity scale iscale={iscale}')

Expand Down Expand Up @@ -339,7 +338,7 @@
"""
output_models = []
for group_id, indices in input_models.group_indices.items():

output_model = self.resample_group(input_models, indices)

if not self.in_memory:
Expand All @@ -357,7 +356,7 @@
# build ModelLibrary as an association from the output files
# this saves memory if there are multiple groups
asn = asn_from_list(output_models, product_name='outlier_i2d')
asn_dict = json.loads(asn.dump()[1]) # serializes the asn and converts to dict
asn_dict = json.loads(asn.dump()[1]) # serializes the asn and converts to dict

Check warning on line 359 in jwst/resample/resample.py

View check run for this annotation

Codecov / codecov/patch

jwst/resample/resample.py#L359

Added line #L359 was not covered by tests
return ModelLibrary(asn_dict, on_disk=True)
# otherwise just build it as a list of in-memory models
return ModelLibrary(output_models, on_disk=False)
Expand Down Expand Up @@ -397,9 +396,11 @@
log.debug(f'Using intensity scale iscale={iscale}')
img.meta.iscale = iscale

inwht = resample_utils.build_driz_weight(img,
weight_type=self.weight_type,
good_bits=self.good_bits)
inwht = resample_utils.build_driz_weight(
img,
weight_type=self.weight_type,
good_bits=self.good_bits,
)
# apply sky subtraction
blevel = img.meta.background.level
if not img.meta.background.subtracted and blevel is not None:
Expand Down Expand Up @@ -446,14 +447,13 @@

return ModelLibrary([output_model,], on_disk=False)


def resample_variance_arrays(self, output_model, input_models):
"""Resample variance arrays from input_models to the output_model.

Variance images from each input model are resampled individually and
added to a weighted sum. If weight_type is 'ivm', the inverse of the
added to a weighted sum. If ``weight_type`` is 'ivm', the inverse of the
resampled read noise variance is used as the weight for all the variance
components. If weight_type is 'exptime', the exposure time is used.
components. If ``weight_type`` is 'exptime', the exposure time is used.

The output_model is modified in place.
"""
Expand Down Expand Up @@ -493,8 +493,10 @@
if rn_var is not None:
mask = (rn_var >= 0) & np.isfinite(rn_var) & (weight > 0)
weighted_rn_var[mask] = np.nansum(
[weighted_rn_var[mask],
rn_var[mask] * weight[mask] * weight[mask]],
[
weighted_rn_var[mask],
rn_var[mask] * weight[mask] * weight[mask]
],
axis=0
)
total_weight_rn_var[mask] += weight[mask]
Expand All @@ -506,8 +508,10 @@
if pn_var is not None:
mask = (pn_var >= 0) & np.isfinite(pn_var) & (weight > 0)
weighted_pn_var[mask] = np.nansum(
[weighted_pn_var[mask],
pn_var[mask] * weight[mask] * weight[mask]],
[
weighted_pn_var[mask],
pn_var[mask] * weight[mask] * weight[mask]
],
axis=0
)
total_weight_pn_var[mask] += weight[mask]
Expand All @@ -517,12 +521,14 @@
if flat_var is not None:
mask = (flat_var >= 0) & np.isfinite(flat_var) & (weight > 0)
weighted_flat_var[mask] = np.nansum(
[weighted_flat_var[mask],
flat_var[mask] * weight[mask] * weight[mask]],
[
weighted_flat_var[mask],
flat_var[mask] * weight[mask] * weight[mask]
],
axis=0
)
total_weight_flat_var[mask] += weight[mask]

del model.meta.iscale
del weight
input_models.shelve(model, i)
Expand All @@ -549,7 +555,6 @@
del weighted_rn_var, weighted_pn_var, weighted_flat_var
del total_weight_rn_var, total_weight_pn_var, total_weight_flat_var


def _resample_one_variance_array(self, name, input_model, output_model):
"""Resample one variance image from an input model.

Expand Down Expand Up @@ -753,7 +758,7 @@
"""

# Insure that the fillval parameter gets properly interpreted for use with tdriz
if util.is_blank(str(fillval)):
if str(fillval).strip() == '':
fillval = 'NAN'
else:
fillval = str(fillval)
Expand Down
16 changes: 5 additions & 11 deletions jwst/resample/resample_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@

import numpy as np
from astropy import units as u
import gwcs
from drizzle.utils import calc_pixmap

from stdatamodels.dqflags import interpret_bit_flags
from stdatamodels.jwst.datamodels.dqflags import pixel
Expand Down Expand Up @@ -124,17 +124,11 @@ def shape_from_bounding_box(bounding_box):
def calc_gwcs_pixmap(in_wcs, out_wcs, shape=None):
""" Return a pixel grid map from input frame to output frame.
"""
if shape:
bb = wcs_bbox_from_shape(shape)
log.debug("Bounding box from data shape: {}".format(bb))
emolter marked this conversation as resolved.
Show resolved Hide resolved
else:
bb = in_wcs.bounding_box
log.debug("Bounding box from WCS: {}".format(in_wcs.bounding_box))

grid = gwcs.wcstools.grid_from_bounding_box(bb)
pixmap = np.dstack(reproject(in_wcs, out_wcs)(grid[0], grid[1]))
if shape is not None and not np.array_equiv(shape, in_wcs.array_shape):
in_wcs = deepcopy(in_wcs)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Is a deepcopy necessary here?

Copy link
Member Author

Choose a reason for hiding this comment

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

I am not sure. It felt to be a safer way because this method sets array_shape property of the WCS object. If there is no drawback in modifying inputs, deepcopy can be omitted.

in_wcs.array_shape = shape

return pixmap
return calc_pixmap(wcs_from=in_wcs, wcs_to=out_wcs)


def reproject(wcs1, wcs2):
Expand Down
2 changes: 1 addition & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ dependencies = [
"astropy>=5.3",
"BayesicFitting>=3.0.1",
"crds>=11.17.14",
"drizzle>=1.15.0",
"drizzle @ git+https://github.com/mcara/drizzle.git@redesign-api",
Copy link
Collaborator

Choose a reason for hiding this comment

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

remember to change to @main before merge

"gwcs>=0.21.0,<0.23.0",
"numpy>=1.22,<2.0",
"opencv-python-headless>=4.6.0.66",
Expand Down
Loading