From 34ccf65a1cc594c30bffb9521a0d81fee807e781 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Tue, 23 Jan 2024 12:42:24 -0500 Subject: [PATCH] Make code run --- drizzle/resample.py | 74 +++++++++++++++++++++++++--------------- drizzle/utils.py | 83 ++++++++++++++++++++++++++++++++------------- 2 files changed, 106 insertions(+), 51 deletions(-) diff --git a/drizzle/resample.py b/drizzle/resample.py index db276b3c..18f42580 100644 --- a/drizzle/resample.py +++ b/drizzle/resample.py @@ -42,10 +42,10 @@ def __init__(self, n_images=1, kernel="square", "lanczos2", and "lanczos3". The square kernel is the default. out_shape : tuple, None, optional - Shape of the output images (context image will have a third - dimension of size proportional to the number of input images). - This parameter is helpful when neither ``out_sci``, ``out_wht``, - nor ``out_ctx`` images are provided. + Shape (`numpy` order ``(Ny, Nx)``) of the output images (context + image will have a third dimension of size proportional to the number + of input images). This parameter is helpful when neither + ``out_sci``, ``out_wht``, nor ``out_ctx`` images are provided. fillval: float, None, str, optional The value of output pixels that did not have contributions from @@ -53,8 +53,8 @@ def __init__(self, n_images=1, kernel="square", ``"INDEF"`` and ``out_sci`` is provided, the values of ``out_sci`` will not be modified. When ``fillval`` is either `None` or ``"INDEF"`` and ``out_sci`` is **not provided**, the values of - ``out_sci`` will be initialized to 0. If ``fillval`` is a number - of a string that can be converted to a number, then the output + ``out_sci`` will be initialized to `numpy.nan`. If ``fillval`` + is a string that can be converted to a number, then the output pixels with no contributions from input images will be set to this ``fillval`` value. @@ -74,8 +74,6 @@ def __init__(self, n_images=1, kernel="square", Subsequent calls hold intermediate results. """ - if not isinstance(ctx_id, numbers.Integral) or ctx_id < 0: - raise TypeError("'ctx_id' must be a non-negative integer number") self.ctx_id = -1 # the ID of the *last* image to be resampled self.exptime = 0.0 @@ -85,6 +83,7 @@ def __init__(self, n_images=1, kernel="square", if fillval is None: fillval = "INDEF" + elif isinstance(str, fillval): fillval = fillval.strip() if fillval == "": @@ -94,6 +93,10 @@ def __init__(self, n_images=1, kernel="square", else: float(fillval) fillval = str(fillval) + + if out_sci is None and fillval == "INDEF": + fillval = "NaN" + self.fillval = fillval if (out_shape is None and out_sci is None and out_wht is None and @@ -104,18 +107,16 @@ def __init__(self, n_images=1, kernel="square", shapes = set() - if out_sci is None: - if fillval == "INDEF": - raise ValueError( - "Fill value 'INDEF' is not allowed when 'out_sci' is None." - ) - else: + if out_sci is not None: + out_sci = np.asarray(out_sci, dtype=np.float32) shapes.add(out_sci.shape) if out_wht is not None: + out_wht = np.asarray(out_wht, dtype=np.float32) shapes.add(out_wht.shape) if out_ctx is not None: + out_ctx = np.asarray(out_ctx, dtype=np.int32) if len(out_ctx.shape) == 3: shapes.add(out_ctx.shape[1:]) else: @@ -144,23 +145,30 @@ def __init__(self, n_images=1, kernel="square", "provided." ) - self._alloc_output_arrays(n_images=n_images) + self._alloc_output_arrays( + out_shape=self.out_shape, + n_images=n_images, + out_sci=out_sci, + out_wht=out_wht, + out_ctx=out_ctx, + ) - def _alloc_output_arrays(self, n_images): + def _alloc_output_arrays(self, out_shape, n_images, out_sci, out_wht, + out_ctx): # allocate arrays as needed: if out_sci is None: - self.out_sci = np.empty(self.out_shape, dtype=np.float32) - self.out_sci.fill(float(self.fillval)) + self.out_sci = np.empty(out_shape, dtype=np.float32) + self.out_sci.fill(self.fillval) if out_wht is None: - self.out_wht = np.ones(self.out_shape, dtype=np.float32) + self.out_wht = np.zeros(out_shape, dtype=np.float32) if out_ctx is None: n_ctx_planes = (n_images - 1) // CTX_PLANE_BITS + 1 if n_ctx_planes == 1: - ctx_shape = self.out_shape + ctx_shape = out_shape else: - ctx_shape = (n_ctx_planes, ) + self.out_shape + ctx_shape = (n_ctx_planes, ) + out_shape self.out_ctx = np.zeros(ctx_shape, dtype=np.int32) def _increment_ctx_id(self): @@ -268,6 +276,12 @@ def add_image(self, data, exptime, pixmap, scale=1.0, self.exptime += exptime data = np.asarray(data, dtype=np.float32) + pixmap = np.asarray(pixmap, dtype=np.float64) + + if pixmap.shape[:2] != data.shape: + raise ValueError( + "'pixmap' shape is not consistent with 'data' shape." + ) # TODO: this is code that would enable initializer to not need output # image shape at all and set output image shape based on @@ -292,14 +306,16 @@ def add_image(self, data, exptime, pixmap, scale=1.0, if ymin is None or ymin < 0: ymin = 0 - if xmax > self.out_shape[1] - 1: + if xmax is None or xmax > self.out_shape[1] - 1: xmax = self.out_shape[1] - 1 - if ymax > self.out_shape[0] - 1: + if ymax is None or ymax > self.out_shape[0] - 1: ymax = self.out_shape[0] - 1 if weight_map is not None: weight_map = np.asarray(weight_map, dtype=np.float32) + else: # TODO: this should not be needed after C code modifications + weight_map = np.ones_like(data) pixmap = np.asarray(pixmap, dtype=np.float64) @@ -309,14 +325,20 @@ def add_image(self, data, exptime, pixmap, scale=1.0, # TODO: While drizzle code in cdrizzlebox.c supports weight_map=None, # cdrizzleapi.c does not. It should be modified to support this # for performance reasons. + if len(self.out_ctx.shape) == 2: + assert plane_no == 0 + ctx_plane = self.out_ctx + else: + ctx_plane = self.out_ctx[plane_no] + _vers, nmiss, nskip = cdrizzle.tdriz( data, weight_map, pixmap, self.out_sci, self.out_wht, - self.out_ctx[plane_no], - uniqid=id_in_plane, + ctx_plane, + uniqid=id_in_plane + 1, xmin=xmin, xmax=xmax, ymin=ymin, @@ -334,8 +356,6 @@ def add_image(self, data, exptime, pixmap, scale=1.0, return nmiss, nskip -# pixmap = calc_pixmap.calc_pixmap(blot_wcs, source_wcs) -# pix_ratio = source_wcs.pscale / blot_wcs.pscale def blot_image(data, pixmap, pix_ratio, exptime, output_pixel_shape, interp='poly5', sinscl=1.0): """ diff --git a/drizzle/utils.py b/drizzle/utils.py index 911c5464..4748689d 100644 --- a/drizzle/utils.py +++ b/drizzle/utils.py @@ -1,52 +1,87 @@ import numpy as np -def calc_pixmap(first_wcs, second_wcs): +def calc_pixmap(wcs_from, wcs_to, estimate_pixel_scale_ratio=False, + refpix_from=None, refpix_to=None): """ - Calculate a mapping between the pixels of two images. + Calculate a mapping between the pixels of two images. Pixel scale ratio, + when requested, is computed near the centers of the bounding box + (a property of the WCS object) or near ``refpix_*`` coordinates + if supplied. Parameters ---------- - first_wcs : wcs + wcs_from : wcs A WCS object representing the coordinate system you are - converting from + converting from. This object *must* have ``pixel_shape`` property + defined. - seond_wcs : wcs + wcs_to : wcs A WCS object representing the coordinate system you are - converting to + converting to. + + estimate_pixel_scale_ratio : bool, optional + Estimate the ratio of "to" to "from" WCS pixel scales. + + refpix_from : numpy.ndarray, tuple, list + Image coordinates of the reference pixel near which pixel scale should + be computed in the "from" image. In FITS WCS this could be, for example, + the value of CRPIX of the ``wcs_from`` WCS. + + refpix_to : numpy.ndarray, tuple, list + Image coordinates of the reference pixel near which pixel scale should + be computed in the "to" image. In FITS WCS this could be, for example, + the value of CRPIX of the ``wcs_to`` WCS. Returns ------- - A three dimensional array representing the transformation between - the two. The last dimension is of length two and contains the x and - y coordinates of a pixel center, repectively. The other two coordinates - correspond to the two coordinates of the image the first WCS is from. - """ + pixmap : numpy.ndarray + A three dimensional array representing the transformation between + the two. The last dimension is of length two and contains the x and + y coordinates of a pixel center, repectively. The other two coordinates + correspond to the two coordinates of the image the first WCS is from. - first_naxis1, first_naxis2 = first_wcs.pixel_shape + pixel_scale_ratio : float + Estimate the ratio of "to" to "from" WCS pixel scales. This value is + returned only when ``estimate_pixel_scale_ratio`` is `True`. + """ # We add one to the pixel co-ordinates before the transformation and subtract # it afterwards because wcs co-ordinates are one based, while pixel co-ordinates # are zero based, The result is the final values in pixmap give a tranformation # between the pixel co-ordinates in the first image to pixel co-ordinates in the # co-ordinate system of the second. - one = np.ones(2, dtype='float64') - idxmap = np.indices((first_naxis1, first_naxis2), dtype='float64') - idxmap = idxmap.transpose() + one + if wcs_from.pixel_shape is None: + raise ValueError('The "from" WCS must have pixel_shape property set.') + y, x = np.indices(wcs_from.pixel_shape, dtype=np.float64) + x, y = wcs_to.invert(*wcs_from(x, y)) + pixmap = np.dstack([x, y]) + if estimate_pixel_scale_ratio: + pscale_ratio = (_estimate_pixel_scale(wcs_to, refpix_to) / + _estimate_pixel_scale(wcs_from, refpix_from)) + return pixmap, pscale_ratio - idxmap = idxmap.reshape(first_naxis2 * first_naxis1, 2) + return pixmap - worldmap = first_wcs.all_pix2world(idxmap, 1) - if second_wcs.sip is None: - pixmap = second_wcs.wcs_world2pix(worldmap, 1) +def _estimate_pixel_scale(wcs, refpix): + # estimate pixel scale (in rad) using approximate algorithm + # from https://trs.jpl.nasa.gov/handle/2014/40409 + if refpix is None: + if wcs.bounding_box is None: + refpix = np.ones(wcs.pixel_n_dim) + else: + refpix = np.mean(wcs.bounding_box, axis=-1) else: - pixmap = second_wcs.all_world2pix(worldmap, 1) + refpix = np.asarray(refpix) - pixmap = pixmap.reshape(first_naxis2, first_naxis1, 2) - pixmap = pixmap - one - - return pixmap + l1, phi1 = np.deg2rad(wcs.__call__(*(refpix - 0.5))) + l2, phi2 = np.deg2rad(wcs.__call__(*(refpix + [-0.5, 0.5]))) + l3, phi3 = np.deg2rad(wcs.__call__(*(refpix + 0.5))) + l4, phi4 = np.deg2rad(wcs.__call__(*(refpix + [0.5, -0.5]))) + area = np.abs(0.5 * ((l4 - l2) * (np.sin(phi1) - np.sin(phi3)) + + (l1 - l3) * (np.sin(phi2) - np.sin(phi4)))) + return np.sqrt(area)