Make code run

drizzle/resample.py

@@ -42,19 +42,19 @@ 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
             input images' pixels. When ``fillval`` is either `None` or
             ``"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):
     """

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)