Fix NDCube.rebin for lookup-table extra coords

changelog/982.bugfix.1.rst

@@ -0,0 +1 @@
+Interpolating a `~ndcube.extra_coords.TimeTableCoordinate`, e.g. when rebinning, now keeps sub-microsecond precision and the coordinate's ``reference_time``.

changelog/982.bugfix.2.rst

@@ -0,0 +1 @@
+`~ndcube.extra_coords.QuantityTableCoordinate.interpolate` now accepts grids of different lengths, so `~ndcube.NDCube.rebin` no longer raises "New array grids must all be same shape" for a `~astropy.units.Quantity` lookup table spanning several axes that are rebinned to different lengths.

changelog/982.bugfix.rst

@@ -0,0 +1 @@
+Fixed `~ndcube.NDCube.rebin` for cubes with lookup-table extra coords: the extra coords are now sampled at the centre of each new pixel, matching the rebinned WCS, instead of at the first pixel of each bin, and rebinning no longer raises a ``TypeError`` when every axis of the rebinned cube has the same length (always the case for 1-D cubes).

ndcube/extra_coords/extra_coords.py

@@ -493,12 +493,11 @@ def resample(self, factor, offset=0, ndcube=None, **kwargs):
             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_coord = coord.interpolate(*[new_grids[axis] for axis in array_axes], **kwargs)
             new_ec.add(coord.names, array_axes, new_coord, physical_types=coord.physical_types)
         return new_ec
 

ndcube/extra_coords/table_coord.py

@@ -448,9 +448,8 @@ def interpolate(self, *new_array_grids, **kwargs):
         new_array_grids: array-like
             The array index values at which the the new values of the coords
             are desired. An array grid must be provided as a separate arg
-            for each array dimension and corresponding elements in all arrays
-            represent a single location in the pixel grid. Therefore, array grids
-            must all have the same shape.
+            for each array dimension. Each table is interpolated independently
+            along its own grid, so the grids can have different lengths.
 
         kwargs
             All remaining kwargs are passed to underlying interpolation function.
@@ -468,8 +467,6 @@ def interpolate(self, *new_array_grids, **kwargs):
         if len(new_array_grids) != ndim:
             raise ValueError(
                 f"A new array grid must be given for each array axis/table, i.e. {ndim}")
-        if any(new_grid.shape != new_array_grids[0].shape for new_grid in new_array_grids):
-            raise ValueError("New array grids must all be same shape.")
         # Build array grids for non-interpolated table.
         old_array_grids = tuple(np.arange(d) for d in self.shape)
         # Iterate through tables and interpolate each.
@@ -796,12 +793,13 @@ def interpolate(self, new_array_grids, **kwargs):
             raise ValueError("Cannot interpolate a scalar TimeTableCoordinate.")
         # Build pixel grids for current TimeTableCoord.
         old_array_grids = np.arange(len(self.table))
-        # Interpolate using MJD format and convert back to a Time object.
-        new_table = np.interp(new_array_grids, old_array_grids, self.table.mjd, **kwargs)
-        new_table = Time(new_table, scale=self.table.scale, format="mjd")
-        new_table.format = self.table.format
+        # Interpolate offsets from the first time; absolute MJD floats lose sub-microsecond precision.
+        origin = self.table[0]
+        deltas = (self.table - origin).to_value(u.s)
+        new_table = origin + np.interp(new_array_grids, old_array_grids, deltas, **kwargs) * u.s
         # Rebuild new TimeTableCoord and return.
-        new_coord = type(self)(new_table, names=self.names, physical_types=self.physical_types)
+        new_coord = type(self)(new_table, names=self.names, physical_types=self.physical_types,
+                               reference_time=self.reference_time)
         new_coord._dropped_world_dimensions = self._dropped_world_dimensions
         return new_coord
 

ndcube/extra_coords/tests/test_lookup_table_coord.py

@@ -658,6 +658,15 @@ def test_time_interpolate(lut_1d_time):
     assert_lutc_ancilliary_data_same(output, lutc)
 
 
+def test_time_interpolate_keeps_precision_and_reference_time():
+    reference = Time("2020-01-01T00:00:00", scale="utc")
+    times = Time("2026-01-01T00:00:00", scale="utc") + np.arange(4) * 10 * u.ns
+    coord = TimeTableCoordinate(times, reference_time=reference)
+    new = coord.interpolate(np.arange(4.))
+    assert new.reference_time == reference
+    assert u.allclose((new.table - times).to(u.ns), 0 * u.ns, atol=1 * u.ns)
+
+
 def test_skycoord_interpolate_no_mesh(lut_2d_skycoord_no_mesh):
     lutc = lut_2d_skycoord_no_mesh
     new_array_grids = np.meshgrid(np.arange(0.5, 2), np.arange(0, 3))
@@ -696,10 +705,6 @@ def test_quantity_interpolate_errors():
         qtc.interpolate(np.ones(1))
     assert "A new array grid must be given for each array axis" in str(ei)
 
-    with pytest.raises(ValueError) as ei:
-        qtc.interpolate(np.ones(1), np.ones(2))
-    assert "New array grids must all be same shape." in str(ei)
-
 
 def test_skycoord_interpolate_error(lut_2d_skycoord_mesh):
     sctc = SkyCoordTableCoordinate(SkyCoord(1 * u.deg, 2 * u.deg))

ndcube/ndcube.py

@@ -1392,7 +1392,9 @@ def my_propagate(uncertainty, data, mask, **kwargs):
         new_cube._global_coords = self._global_coords
         # Reconstitute extra coords
         if not self.extra_coords.is_empty:
-            new_cube._extra_coords = self.extra_coords.resample(bin_shape, ndcube=new_cube)
+            # Lookup tables are sampled at bin centres; WCS wrappers take pixel-edge offsets.
+            offset = (bin_shape - 1) / 2 if self.extra_coords._lookup_tables else 0
+            new_cube._extra_coords = self.extra_coords.resample(bin_shape, offset=offset, ndcube=new_cube)
 
         return new_cube
 

ndcube/tests/test_ndcube_reproject_and_rebin.py

@@ -157,7 +157,7 @@ def test_rebin(ndcube_3d_l_ln_lt_ectime, bin_shape):
     assert u.allclose(output_spec, expected_spec)
     assert output_time.scale == expected_time.scale
     assert output_time.format == expected_time.format
-    assert np.allclose(output_time.mjd, expected_time.mjd)
+    assert np.allclose(output_time.mjd, expected_time.mjd, rtol=0, atol=1e-8)
 
 
 def test_rebin_dask(ndcube_2d_dask):
@@ -325,3 +325,18 @@ def test_rebin_specutils():
     spec = Spectrum(flux=y, spectral_axis=x, bin_specification='centers', mask=x > 2000*u.nm)
     output = spec.rebin((10,), operation=np.sum, operation_ignores_mask=False)
     assert output.shape == (400,)
+
+
+def test_rebin_1d_extra_coords(ndcube_1d_l):
+    ndcube_1d_l.extra_coords.add("wave", 0, ndcube_1d_l.axis_world_coords_values()[0])
+    output = ndcube_1d_l.rebin((2,))
+    assert u.allclose(output.axis_world_coords_values(wcs=output.extra_coords)[0],
+                      output.axis_world_coords_values()[0])
+
+
+def test_rebin_mesh_extra_coords(ndcube_2d_ln_lt):
+    ndcube_2d_ln_lt.extra_coords.add(("a", "b"), (0, 1), (np.arange(10) * u.m, np.arange(12) * u.m))
+    output = ndcube_2d_ln_lt.rebin((5, 4))
+    a, b = output.axis_world_coords(wcs=output.extra_coords)
+    assert u.allclose(a, [2, 7] * u.m)
+    assert u.allclose(b, [1.5, 5.5, 9.5] * u.m)