Implement Superpixel on NDCube

changelog/450.feature.1.rst

@@ -0,0 +1 @@
+Implement a new `ndcube.NDCube.rebin` method to combine integer numbers of pixels along each axis into larger pixels.

changelog/450.feature.2.rst

@@ -0,0 +1 @@
+Implement new property, `ndcube.ExtraCoords.is_empty` that returns ``True`` if the object has not extra coords.  Otherwise return ``False``.

changelog/450.feature.3.rst

@@ -0,0 +1 @@
+Add new methods to interpolate lookup table coordinates: `ndcube.extra_coords.table_coord.QuantityTableCoordinate.interpolate`, `ndcube.extra_coords.table_coord.SkyCoordTableCoordinate.interpolate`, `ndcube.extra_coords.table_coord.TimeTableCoordinate.interpolate`, `ndcube.extra_coords.table_coord.MultipleTableCoordinate.interpolate`

changelog/450.feature.4.rst

@@ -0,0 +1 @@
+Add `ndcube.ExtraCoords.resample` method to resample extra coordinates by a certain factor in each array dimension.

ndcube/conftest.py

@@ -4,6 +4,7 @@
 """
 import logging
 
+import astropy.nddata
 import astropy.units as u
 import dask.array
 import numpy as np
@@ -58,11 +59,15 @@ def gen_ndcube_3d_l_ln_lt_ectime(wcs_3d_lt_ln_l, time_axis, time_base, global_co
     wcs_3d_lt_ln_l.array_shape = shape
     data_cube = data_nd(shape)
     mask = data_cube < 0
+    meta = {"message": "hello world"}
+    unit = u.ph
     extra_coords = time_extra_coords(shape, time_axis, time_base)
     cube = NDCube(data_cube,
                   wcs_3d_lt_ln_l,
                   mask=mask,
-                  uncertainty=data_cube)
+                  uncertainty=data_cube,
+                  meta=meta,
+                  unit=unit)
     cube._extra_coords = extra_coords
 
     if global_coords:
@@ -468,6 +473,41 @@ def ndcube_2d_ln_lt(wcs_2d_lt_ln):
     return NDCube(data_cube, wcs=wcs_2d_lt_ln)
 
 
+@pytest.fixture
+def ndcube_2d_ln_lt_uncert(wcs_2d_lt_ln):
+    shape = (10, 12)
+    data_cube = data_nd(shape)
+    uncertainty = astropy.nddata.StdDevUncertainty(data_cube * 0.1)
+    cube = NDCube(data_cube, wcs=wcs_2d_lt_ln, uncertainty=uncertainty)
+    return cube
+
+
+@pytest.fixture
+def ndcube_2d_ln_lt_mask_uncert(wcs_2d_lt_ln):
+    shape = (10, 12)
+    data_cube = data_nd(shape)
+    uncertainty = astropy.nddata.StdDevUncertainty(data_cube * 0.1)
+    mask = np.zeros(shape, dtype=bool)
+    mask[1, 1] = True
+    mask[2, 0] = True
+    mask[3, 3] = True
+    mask[4:6, :4] = True
+    cube = NDCube(data_cube, wcs=wcs_2d_lt_ln, uncertainty=uncertainty, mask=mask)
+    return cube
+
+
+@pytest.fixture
+def ndcube_2d_ln_lt_uncert_ec(wcs_2d_lt_ln):
+    shape = (4, 9)
+    data_cube = data_nd(shape)
+    uncertainty = astropy.nddata.StdDevUncertainty(data_cube * 0.1)
+    cube = NDCube(data_cube, wcs=wcs_2d_lt_ln, uncertainty=uncertainty)
+    cube.extra_coords.add(
+        "time", 0,
+        Time("2000-01-01 00:00", scale="utc") + Timedelta(np.arange(shape[0])*60, format="sec"))
+    return cube
+
+
 @pytest.fixture
 def ndcube_2d_ln_lt_units(wcs_2d_lt_ln):
     shape = (10, 12)
@@ -480,11 +520,11 @@ def ndcube_2d_dask(wcs_2d_lt_ln):
     shape = (8, 4)
     chunks = 2
     data = data_nd(shape).astype(float)
-    da = dask.array.from_array(data, chunks=chunks)
+    da = dask.array.asarray(data, chunks=chunks)
     mask = np.zeros(shape, dtype=bool)
-    da_mask = dask.array.from_array(mask, chunks=chunks)
+    da_mask = dask.array.asarray(mask, chunks=chunks)
     uncert = data * 0.1
-    da_uncert = StdDevUncertainty(dask.array.from_array(uncert, chunks=chunks))
+    da_uncert = StdDevUncertainty(dask.array.asarray(uncert, chunks=chunks))
     return NDCube(da, wcs=wcs_2d_lt_ln, uncertainty=da_uncert, mask=da_mask, unit=u.J)
 
 

ndcube/extra_coords/extra_coords.py

@@ -9,10 +9,11 @@
 from astropy.time import Time
 from astropy.wcs import WCS
 from astropy.wcs.wcsapi import BaseHighLevelWCS
+from astropy.wcs.wcsapi.high_level_wcs_wrapper import HighLevelWCSWrapper
 from astropy.wcs.wcsapi.wrappers.sliced_wcs import SlicedLowLevelWCS, sanitize_slices
 
 from ndcube.utils.wcs import convert_between_array_and_pixel_axes
-from ndcube.wcs.wrappers import CompoundLowLevelWCS
+from ndcube.wcs.wrappers import CompoundLowLevelWCS, ResampledLowLevelWCS
 
 from .table_coord import (BaseTableCoordinate, MultipleTableCoordinate, QuantityTableCoordinate,
                           SkyCoordTableCoordinate, TimeTableCoordinate)
@@ -53,7 +54,7 @@ def add(self,
         name : `str` or sequence of `str`
             The name(s) for these world coordinate(s).
         array_dimension : `int` or `tuple` of `int`
-            The pixel dimension(s), in the array, to which this lookup table corresponds.
+            The array dimension(s), in the array, to which this lookup table corresponds.
         lookup_table : `object` or sequence of `object`
             The lookup table. A `BaseTableCoordinate <.table_coord>` subclass or anything
             that can instantiate one, i.e. currently a `~astropy.time.Time`,
@@ -96,6 +97,11 @@ def wcs(self) -> BaseHighLevelWCS:
 
         """
 
+    @property
+    @abc.abstractproperty
+    def is_empty(self):
+        """Return True if no extra coords present, else return False."""
+
     @abc.abstractmethod
     def __getitem__(self, item: Union[str, int, slice, Iterable[Union[str, int, slice]]]) -> "ExtraCoordsABC":
         """
@@ -218,7 +224,7 @@ def add(self, name, array_dimension, lookup_table, physical_types=None, **kwargs
 
         # Sort the LUTs so that the mapping and the wcs are ordered in pixel dim order
         self._lookup_tables = list(sorted(self._lookup_tables,
-                                          key=lambda x: x[0] if isinstance(x[0], int) else x[0][0]))
+                                          key=lambda x: x[0] if isinstance(x[0], Integral) else x[0][0]))
 
     @property
     def _name_lut_map(self):
@@ -305,6 +311,14 @@ def wcs(self, wcs):
 
         self._wcs = wcs
 
+    @property
+    def is_empty(self):
+        # docstring in ABC
+        if not self._wcs and not self._lookup_tables:
+            return True
+        else:
+            return False
+
     def _getitem_string(self, item):
         """
         Slice the Extracoords based on axis names.
@@ -406,6 +420,72 @@ def dropped_world_dimensions(self):
 
         return dict()
 
+    def resample(self, factor, offset=0, ndcube=None, **kwargs):
+        """
+        Resample all extra coords by given factors in array-index-space.
+
+        One resample factor must be supplied for each array axis in array-axis order.
+        Kwargs are passed to `numpy.interp`.
+
+        Parameters
+        ----------
+        factor: `ìnt`, `float`, or iterable thereof.
+            The factor by which each array axis is resampled.
+            If scalar, same factor is applied to all axes.
+            Otherwise a factor for each axis must be provided.
+
+        offset: `int` `float` of iterable therefore.
+            The location on the underlying grid which corresponds
+            to the zeroth element after resampling. If iterable, must have an entry
+            for each dimension.  If a scalar, the grid will be
+            shifted by the same amount in all dimensions.
+
+        ndcube: `~ndcube.NDCube`
+            The NDCube instance with which the output ExtraCoords object is associated.
+
+        Returns
+        -------
+        new_ec: `~ndcube.extra_coords.ExtraCoords`
+            A new ExtraCoords object holding the interpolated coords.
+        """
+        new_ec = type(self)(ndcube)
+        if self.is_empty:
+            return new_ec
+        cube_shape = self._ndcube.data.shape
+        ndim = len(cube_shape)
+        if np.isscalar(factor):
+            factor = [factor] * ndim
+        if len(factor) != ndim:
+            raise ValueError(
+                "factor must be scalar or an iterable with length equal to number of cube "
+                f"dimensions: len(factor) = {len(factor)}; No. cube dimensions = {ndim}.")
+        if np.isscalar(offset):
+            offset = [offset] * ndim
+        if len(offset) != ndim:
+            raise ValueError(
+                "offset must be scalar or an iterable with length equal to number of cube "
+                f"dimensions: len(offset) = {len(offset)}; No. cube dimensions = {ndim}.")
+        # If ExtraCoords object built on WCS, resample using WCS insfrastructure
+        if self._wcs is not None:
+            new_ec.wcs = HighLevelWCSWrapper(ResampledLowLevelWCS(self._wcs.low_level_wcs,
+                                                                  factor, offset))
+            return new_ec
+        # Else interpolate the lookup table coordinates.
+        factor = np.asarray(factor)
+        new_grids = []
+        for c, d, f in zip(offset, cube_shape, factor):
+            x = np.arange(c, d+f, f)
+            x = x[x <= d-1]
+            new_grids.append(x)
+        new_grids = np.array(new_grids, dtype=object)
+        for array_axes, coord in self._lookup_tables:
+            if np.isscalar(array_axes):
+                new_coord = coord.interpolate(new_grids[array_axes], **kwargs)
+            else:
+                new_coord = coord.interpolate(*new_grids[np.asarray(array_axes)], **kwargs)
+            new_ec.add(coord.names, array_axes, new_coord, physical_types=coord.physical_types)
+        return new_ec
+
     @property
     def cube_wcs(self):
         """Produce a WCS that describes the associated NDCube with just the extra coords.