From c072859d722cb790a2f627fbeae4f4408025ee1e Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Thu, 14 Sep 2023 04:55:56 -0400 Subject: [PATCH] Fix computation of the photometric keywords PIXAR_* and scale resampled intensities according to actual pixel scale ratio. --- CHANGES.rst | 19 ++- docs/jwst/resample/arguments.rst | 18 ++- jwst/model_blender/blendmeta.py | 8 +- jwst/resample/gwcs_drizzle.py | 17 ++- jwst/resample/resample.py | 198 +++++++++++++++++++++++++++++-- jwst/resample/resample_step.py | 23 ++-- 6 files changed, 255 insertions(+), 28 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 061eaf313d6..b1462cfa228 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -65,7 +65,7 @@ charge_migration - Step was renamed from undersampling_migration. Changed default signal threshold, added efficient routine to flag neighborhood pixels, added new unit test, improved earlier unit tests, updated docs. [#7825] - + cube_build ---------- @@ -141,7 +141,7 @@ pathloss - Fix interpolation error for point source corrections. [#7799] - Update the MIRI LRS fixed-slit correction to default to the center of the slit - when the computed target location is outside the slit. Add the parameter + when the computed target location is outside the slit. Add the parameter "user_slit_loc" to allow specifying the source location to use. [#7806] photom @@ -179,12 +179,23 @@ resample - Update the following exposure time keywords: XPOSURE (EFFEXPTM), DURATION and TELAPSE. [#7793] +- Recognize additional keys in ASDF files that provide ``output_wcs`` for the + resample step. [#7894] + +- Set output image size when ``output_wcs`` is provided based on the largest + coordinates in the bounding box of the ``output_wcs``. [#7894] + +- Completely re-designed computation of the pixel area keywords + ``PIXAR_SR`` and ``PIXAR_A2`` for the resampled image. This change also + results in modified values in the resampled images. New computations + significantly reduce photometric errors. [#7894] + residual_fringe --------------- - Use scipy.interpolate.BSpline instead of astropy.modeling.Spline1D in residual_fringe fitting utils [#7764] - + undersampling_correction ------------------------ @@ -193,7 +204,7 @@ undersampling_correction - Removed directories for undersampling_correction step, as the step has been renamed charge_migration. [#7850] - + 1.11.4 (2023-08-14) =================== diff --git a/docs/jwst/resample/arguments.rst b/docs/jwst/resample/arguments.rst index fb3d1777937..ac597eacae4 100644 --- a/docs/jwst/resample/arguments.rst +++ b/docs/jwst/resample/arguments.rst @@ -62,7 +62,23 @@ image. under the root of the file. The output image size is determined from the bounding box of the WCS (if any). Argument ``output_shape`` overrides computed image size and it is required when output WCS does not have - ``bounding_box`` property set. + ``bounding_box`` property set or if ``pixel_shape`` or ``array_shape`` keys + (see below) are not provided. + + Additional information may be stored under + other keys under the root of the file. Currently, the following keys are + recognized: + + - ``pixel_area``: Indicates average pixel area of the output WCS in + units of steradians. When provided, this value will be used for updating + photometric quantities ``PIXAR_SR`` and ``PIXAR_A2`` of the output image. + If ``pixel_area`` is not provided, the code will attempt to estimate + this value from the WCS. + + - ``pixel_shape``: dimensions of the output image in the order (nx, ny). + Overrides the value of ``array_shape`` if provided. + + - ``array_shape``: shape of the output image in ``numpy`` order: (ny, nx). .. note:: When ``output_wcs`` is specified, WCS-related arguments such as diff --git a/jwst/model_blender/blendmeta.py b/jwst/model_blender/blendmeta.py index bf151a37028..9d68886140a 100644 --- a/jwst/model_blender/blendmeta.py +++ b/jwst/model_blender/blendmeta.py @@ -25,7 +25,7 @@ # Primary functional interface for the code -def blendmodels(product, inputs=None, output=None, verbose=False): +def blendmodels(product, inputs=None, output=None, ignore=None, verbose=False): """ Run main interface for blending metadata from multiple models. @@ -73,6 +73,10 @@ def blendmodels(product, inputs=None, output=None, verbose=False): If provided, update `meta.filename` in the blended `product` to define what file this model will get written out to. + ignore : list of str, None, optional + A list of string the meta attribute names which, if provided, + will show which attributes should not be blended. + verbose : bool, optional [Default: False] Print out additional messages during processing when specified. @@ -122,6 +126,8 @@ def blendmodels(product, inputs=None, output=None, verbose=False): # Start by identifying elements of the model which need to be ignored ignore_list = _build_schema_ignore_list(newmeta._schema) ignore_list += ['meta.wcs'] # Necessary since meta.wcs is not in schema + if ignore: + ignore_list.extend(ignore) # Now assign values from new_hdrs to output_model.meta flat_new_metadata = newmeta.to_flat_dict() diff --git a/jwst/resample/gwcs_drizzle.py b/jwst/resample/gwcs_drizzle.py index 64a4a136d1c..8337da121b4 100644 --- a/jwst/resample/gwcs_drizzle.py +++ b/jwst/resample/gwcs_drizzle.py @@ -120,7 +120,7 @@ def outcon(self, value): self._product.con = value def add_image(self, insci, inwcs, inwht=None, xmin=0, xmax=0, ymin=0, ymax=0, - expin=1.0, in_units="cps", wt_scl=1.0): + expin=1.0, in_units="cps", wt_scl=1.0, iscale=1.0): """ Combine an input image with the output drizzled image. @@ -185,6 +185,11 @@ def add_image(self, insci, inwcs, inwht=None, xmin=0, xmax=0, ymin=0, ymax=0, initialized with wt_scl set to "exptime" or "expsq", the exposure time will be used to set the weight scaling and the value of this parameter will be ignored. + + iscale : float, optional + A scale factor to be applied to pixel intensities of the + input image before resampling. + """ if self.wt_scl == "exptime": wt_scl = expin @@ -197,7 +202,8 @@ def add_image(self, insci, inwcs, inwht=None, xmin=0, xmax=0, ymin=0, ymax=0, dodrizzle(insci, inwcs, inwht, self.outwcs, self.outsci, self.outwht, self.outcon, expin, in_units, wt_scl, uniqid=self.uniqid, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, - pixfrac=self.pixfrac, kernel=self.kernel, fillval=self.fillval) + iscale=iscale, pixfrac=self.pixfrac, kernel=self.kernel, + fillval=self.fillval) def increment_id(self): """ @@ -227,7 +233,7 @@ def increment_id(self): def dodrizzle(insci, input_wcs, inwht, output_wcs, outsci, outwht, outcon, expin, in_units, wt_scl, uniqid=1, xmin=0, xmax=0, ymin=0, ymax=0, - pixfrac=1.0, kernel='square', fillval="INDEF"): + iscale=1.0, pixfrac=1.0, kernel='square', fillval="INDEF"): """ Low level routine for performing 'drizzle' operation on one image. @@ -303,6 +309,10 @@ def dodrizzle(insci, input_wcs, inwht, output_wcs, outsci, outwht, outcon, no maximum will be set in the y dimension (all pixels in a column of the input image will be resampled). + iscale : float, optional + A scale factor to be applied to pixel intensities of the + input image before resampling. + pixfrac : float, optional The fraction of a pixel that the pixel flux is confined to. The default value of 1 has the pixel flux evenly spread across the image. @@ -389,6 +399,7 @@ def dodrizzle(insci, input_wcs, inwht, output_wcs, outsci, outwht, outcon, uniqid=uniqid, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + scale=iscale, pixfrac=pixfrac, kernel=kernel, in_units=in_units, diff --git a/jwst/resample/resample.py b/jwst/resample/resample.py index 45e739b72bc..a9df5a503b3 100644 --- a/jwst/resample/resample.py +++ b/jwst/resample/resample.py @@ -3,6 +3,7 @@ import numpy as np from drizzle import util from drizzle import cdrizzle +from spherical_geometry.polygon import SphericalPolygon from stdatamodels.jwst import datamodels @@ -76,6 +77,8 @@ def __init__(self, input_models, output=None, single=False, blendheaders=True, self.weight_type = wht_type self.good_bits = good_bits self.in_memory = kwargs.get('in_memory', True) + self.input_pixscale0 = None # computed pixel scale of the first image (deg) + self._recalc_pscale_ratio = pscale is not None log.info(f"Driz parameter kernel: {self.kernel}") log.info(f"Driz parameter pixfrac: {self.pixfrac}") @@ -100,6 +103,15 @@ def __init__(self, input_models, output=None, single=False, blendheaders=True, if output_shape is not None: self.output_wcs.array_shape = output_shape[::-1] + if output_wcs.pixel_area is None: + output_pix_area = _compute_image_pixel_area(self.output_wcs) + if output_pix_area is None: + raise ValueError( + "Unable to compute output pixel area from 'output_wcs'." + ) + else: + output_pix_area = output_wcs.pixel_area + else: # Define output WCS based on all inputs, including a reference WCS: self.output_wcs = resample_utils.make_output_wcs( @@ -113,6 +125,21 @@ def __init__(self, input_models, output=None, single=False, blendheaders=True, crval=crval ) + # Estimate output pixel area in Sr. NOTE: in principle we could + # use the same algorithm as for when output_wcs is provided by the + # user. + tr = self.output_wcs.pipeline[0].transform + output_pix_area = ( + np.deg2rad(tr['cdelt1'].factor.value) * + np.deg2rad(tr['cdelt2'].factor.value) + ) + + if pscale is None: + pscale = np.rad2deg(np.sqrt(output_pix_area)) + log.info(f'Computed output pixel scale: {3600.0 * pscale} arcsec.') + + self.pscale = pscale # in deg + log.debug('Output mosaic size: {}'.format(self.output_wcs.array_shape)) can_allocate, required_memory = datamodels.util.check_memory_allocation( self.output_wcs.array_shape, kwargs['allowed_memory'], datamodels.ImageModel @@ -128,6 +155,10 @@ def __init__(self, input_models, output=None, single=False, blendheaders=True, # update meta data and wcs self.blank_output.update(input_models[0]) self.blank_output.meta.wcs = self.output_wcs + self.blank_output.meta.photometry.pixelarea_steradians = output_pix_area + self.blank_output.meta.photometry.pixelarea_arcsecsq = np.rad2deg( + 3600 * output_pix_area + )**2 self.output_models = ModelContainer(open_models=False) @@ -144,8 +175,18 @@ def blend_output_metadata(self, output_model): # Run fitsblender on output product output_file = output_model.meta.filename + ignore_list = [ + 'meta.photometry.pixelarea_steradians', + 'meta.photometry.pixelarea_arcsecsq', + ] + log.info('Blending metadata for {}'.format(output_file)) - blendmeta.blendmodels(output_model, inputs=self.input_models, output=output_file) + blendmeta.blendmodels( + output_model, + inputs=self.input_models, + output=output_file, + ignore=ignore_list + ) def resample_many_to_many(self): """Resample many inputs to many outputs where outputs have a common frame. @@ -173,6 +214,24 @@ def resample_many_to_many(self): log.info(f"{len(exposure)} exposures to drizzle together") for img in exposure: img = datamodels.open(img) + + input_pixflux_area = img.meta.photometry.pixelarea_steradians + if input_pixflux_area: + img.meta.wcs.array_shape = img.data.shape + input_pixel_area = _compute_image_pixel_area(img.meta.wcs) + if input_pixel_area is None: + raise ValueError( + "Unable to compute input pixel area from WCS of input " + f"image {repr(img.meta.filename)}." + ) + if self.input_pixscale0 is None: + self.input_pixscale0 = np.rad2deg(np.sqrt(input_pixel_area)) + if self._recalc_pscale_ratio: + self.pscale_ratio = self.pscale / self.input_pixscale0 + iscale = np.sqrt(input_pixflux_area / input_pixel_area) + else: + iscale = 1.0 + # TODO: should weight_type=None here? inwht = resample_utils.build_driz_weight( img, @@ -195,6 +254,7 @@ def resample_many_to_many(self): driz.add_image( data, img.meta.wcs, + iscale=iscale, inwht=inwht, xmin=xmin, xmax=xmax, @@ -236,6 +296,23 @@ def resample_many_to_one(self): log.info("Resampling science data") for img in self.input_models: + input_pixflux_area = img.meta.photometry.pixelarea_steradians + if input_pixflux_area: + img.meta.wcs.array_shape = img.data.shape + input_pixel_area = _compute_image_pixel_area(img.meta.wcs) + if input_pixel_area is None: + raise ValueError( + "Unable to compute input pixel area from WCS of input " + f"image {repr(img.meta.filename)}." + ) + if self.input_pixscale0 is None: + self.input_pixscale0 = np.rad2deg(np.sqrt(input_pixel_area)) + if self._recalc_pscale_ratio: + self.pscale_ratio = self.pscale / self.input_pixscale0 + iscale = np.sqrt(input_pixflux_area / input_pixel_area) + else: + iscale = 1.0 + inwht = resample_utils.build_driz_weight(img, weight_type=self.weight_type, good_bits=self.good_bits) @@ -254,6 +331,7 @@ def resample_many_to_one(self): driz.add_image( data, img.meta.wcs, + iscale=iscale, inwht=inwht, xmin=xmin, xmax=xmax, @@ -263,9 +341,9 @@ def resample_many_to_one(self): del data, inwht # Resample variances array in self.input_models to output_model - self.resample_variance_array("var_rnoise", output_model) - self.resample_variance_array("var_poisson", output_model) - self.resample_variance_array("var_flat", output_model) + self.resample_variance_array("var_rnoise", output_model, iscale=iscale) + self.resample_variance_array("var_poisson", output_model, iscale=iscale) + self.resample_variance_array("var_flat", output_model, iscale=iscale) output_model.err = np.sqrt( np.nansum( [ @@ -282,7 +360,7 @@ def resample_many_to_one(self): return self.output_models - def resample_variance_array(self, name, output_model): + def resample_variance_array(self, name, output_model, iscale): """Resample variance arrays from self.input_models to the output_model Resample the ``name`` variance array to the same name in output_model, @@ -335,6 +413,7 @@ def resample_variance_array(self, name, output_model): resampled_variance, outwht, outcon, + iscale=iscale**2, pixfrac=self.pixfrac, kernel=self.kernel, fillval=np.nan, @@ -387,8 +466,8 @@ def update_exposure_times(self, output_model): @staticmethod def drizzle_arrays(insci, inwht, input_wcs, output_wcs, outsci, outwht, outcon, uniqid=1, xmin=0, xmax=0, ymin=0, ymax=0, - pixfrac=1.0, kernel='square', fillval="INDEF", - wtscale=1.0): + iscale=1.0, pixfrac=1.0, kernel='square', + fillval="INDEF", wtscale=1.0): """ Low level routine for performing 'drizzle' operation on one image. @@ -459,6 +538,10 @@ def drizzle_arrays(insci, inwht, input_wcs, output_wcs, outsci, outwht, no maximum will be set in the y dimension (all pixels in a column of the input image will be resampled). + iscale : float, optional + A scale factor to be applied to pixel intensities of the + input image before resampling. + pixfrac : float, optional The fraction of a pixel that the pixel flux is confined to. The default value of 1 has the pixel flux evenly spread across the image. @@ -533,6 +616,7 @@ def drizzle_arrays(insci, inwht, input_wcs, output_wcs, outsci, outwht, uniqid=uniqid, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + scale=iscale, pixfrac=pixfrac, kernel=kernel, in_units="cps", @@ -540,3 +624,103 @@ def drizzle_arrays(insci, inwht, input_wcs, output_wcs, outsci, outwht, wtscale=wtscale, fillstr=fillval ) + + +def _get_boundary_points(xmin, xmax, ymin, ymax, dx=None, dy=None, shrink=0): + """ + xmin, xmax, ymin, ymax - integer coordinates of pixel boundaries + step - distance between points along an edge + shrink - number of pixels by which to reduce `shape` + + Returns a list of points and the area of the rectangle + """ + nx = xmax - xmin + 1 + ny = ymax - ymin + 1 + + if dx is None: + dx = nx + if dy is None: + dy = ny + + if nx - 2 * shrink < 1 or ny - 2 * shrink < 1: + raise ValueError("Image size is too small.") + + sx = max(1, int(np.ceil(nx / dx))) + sy = max(1, int(np.ceil(ny / dy))) + + xmin = shrink + xmax = nx - 1 - shrink + ymin = shrink + ymax = ny - 1 - shrink + + size = 2 * sx + 2 * sy + 1 + x = np.empty(size) + y = np.empty(size) + + x[0:sx] = np.linspace(xmin, xmax, sx, False) + y[0:sx] = ymin + x[sx:sx + sy] = xmax + y[sx:sx + sy] = np.linspace(ymin, ymax, sy, False) + x[sx + sy:2 * sx + sy] = np.linspace(xmax, xmin, sx, False) + y[sx + sy:2 * sx + sy] = ymax + x[2 * sx + sy:2 * sx + 2 * sy] = xmin + y[2 * sx + sy:2 * sx + 2 * sy] = np.linspace(ymax, ymin, sy, False) + x[-1] = xmin + y[-1] = ymin + + area = (xmax - xmin) * (ymax - ymin) + center = (0.5 * (xmin + xmax), 0.5 * (ymin + ymax)) + + return x, y, area, center + + +def _compute_image_pixel_area(wcs): + """ Computes pixel area in steradians. + """ + if wcs.array_shape is None: + raise ValueError("WCS must have array_shape attribute set.") + + valid_polygon = False + spatial_idx = np.where(np.array(wcs.output_frame.axes_type) == 'SPATIAL')[0] + + ny, nx = wcs.array_shape + + ((xmin, xmax), (ymin, ymax)) = wcs.bounding_box + + xmin = max(0, int(xmin + 0.5)) + xmax = min(nx - 1, int(xmax - 0.5)) + ymin = max(0, int(ymin + 0.5)) + ymax = min(ny - 1, int(ymax - 0.5)) + + for shrink in range(30): + try: + x, y, image_area, center = _get_boundary_points( + xmin=xmin, + xmax=xmax, + ymin=ymin, + ymax=ymax, + dx=min(nx // 4, 15), + dy=min(ny // 4, 15), + shrink=shrink + ) + except ValueError: + return None + + world = wcs(x, y) + ra = world[spatial_idx[0]] + dec = world[spatial_idx[1]] + + if (np.all(np.isfinite(ra)) and np.all(np.isfinite(dec))): + valid_polygon = True + break + + if not valid_polygon: + return None + + world = wcs(*center) + wcenter = (world[spatial_idx[0]], world[spatial_idx[1]]) + + sky_area = SphericalPolygon.from_radec(ra, dec, center=wcenter).area() + pix_area = sky_area / image_area + + return pix_area diff --git a/jwst/resample/resample_step.py b/jwst/resample/resample_step.py index 5202c498efa..a3bde264588 100755 --- a/jwst/resample/resample_step.py +++ b/jwst/resample/resample_step.py @@ -140,12 +140,11 @@ def process(self, input): # if pixel_scale exists, it will override pixel_scale_ratio. # calculate the actual value of pixel_scale_ratio based on pixel_scale # because source_catalog uses this value from the header. - if not self.pixel_scale: + if self.pixel_scale is None: model.meta.resample.pixel_scale_ratio = self.pixel_scale_ratio else: - model.meta.resample.pixel_scale_ratio = self.pixel_scale / np.sqrt(model.meta.photometry.pixelarea_arcsecsq) + model.meta.resample.pixel_scale_ratio = resamp.pscale_ratio model.meta.resample.pixfrac = kwargs['pixfrac'] - self.update_phot_keywords(model) if len(result) == 1: result = result[0] @@ -176,17 +175,23 @@ def _load_custom_wcs(asdf_wcs_file, output_shape): with asdf.open(asdf_wcs_file) as af: wcs = deepcopy(af.tree["wcs"]) + wcs.pixel_area = af.tree.get("pixel_area", None) + wcs.array_shape = af.tree.get("pixel_shape", None) + wcs.array_shape = af.tree.get("array_shape", None) - if output_shape is not None or wcs is None: + if output_shape is not None: wcs.array_shape = output_shape[::-1] + wcs.pixel_shape = output_shape elif wcs.pixel_shape is not None: wcs.array_shape = wcs.pixel_shape[::-1] + elif wcs.array_shape is not None: + wcs.pixel_shape = wcs.array_shape[::-1] elif wcs.bounding_box is not None: wcs.array_shape = tuple( - int(axs[1] - axs[0] + 0.5) + int(axs[1] + 0.5) for axs in wcs.bounding_box.bounding_box(order="C") ) - elif wcs.array_shape is None: + else: raise ValueError( "Step argument 'output_shape' is required when custom WCS " "does not have neither of 'array_shape', 'pixel_shape', or " @@ -195,12 +200,6 @@ def _load_custom_wcs(asdf_wcs_file, output_shape): return wcs - def update_phot_keywords(self, model): - """Update pixel scale keywords""" - if model.meta.photometry.pixelarea_steradians is not None: - model.meta.photometry.pixelarea_steradians *= model.meta.resample.pixel_scale_ratio**2 - if model.meta.photometry.pixelarea_arcsecsq is not None: - model.meta.photometry.pixelarea_arcsecsq *= model.meta.resample.pixel_scale_ratio**2 def get_drizpars(self, ref_filename, input_models): """